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ABSTRACT 

The study of how stars distribute themselves around a massive black hole (MBH) in the center of a galaxy is 
an important prerequisite for the understanding of many galactic-center processes. These include the observed 
overabundance of point X-ray sources at the Galactic center, the prediction of rates and characteristics of 
tidal disruptions of extended stars by the MBH and of inspirals of compact stars into the MBH, the latter being 
events of high importance for the future space borne gravitational wave interferometer LISA. In relatively small 
galactic nuclei, hosting MBHs with masses in the range 10 5 -10 7 M Q , the single most important dynamical 
process is 2-body relaxation. It induces the formation of a steep density cusp around the MBH and strong 
mass segregation, as more massive stars lose energy to lighter ones and drift to the central regions. Using 
a spherical stellar dynamical Monte-Carlo code, we simulate the long-term relaxational evolution of galactic 
nucleus models with a spectrum of stellar masses. Our focus is the concentration of stellar black holes to the 
immediate vicinity of the MBH. We quantify this mass segregation for a variety of galactic nucleus models and 
discuss its astrophysical implications. Special attention is given to models developed to match the conditions 
in the Milky Way nucleus; we examine the presence of compact objects in connection to recent high-resolution 
X-ray observations. 

Subject headings: Galaxies: Nuclei, Black Holes, Star Clusters — Methods: A/-Body Simulations, Stellar 
Dynamics — Gravitational Waves 



1. INTRODUCTION 

Massive black holes (MBHs), with masses ranging from a 
few 10 4 M Q to a few 10 9 M© are probably present in the cen- 
ters of most galaxies. The most compelling line of evidence is 
based on measurements of the kinematics of gas and st ars in 
the central regions of nearbv galaxies (e.g., Barth 2004; Kor- 
mendv l2004URichstonel2004tlFerrarese & Fordl2005L for re- 
cent reviews). The inferred masses of the central dark objects 
correlate with different properties of the host galaxy, proba- 
bly most tightly and most fundamentally with the overall ve- 
locity dispersion of the sp heroidal stellar comp onent of the 
galaxy (M-a relation. seelGebhardt et al.ll2000t Ferrar ese & 
Merritt 120001: iTremaine et alJl2002t iNovak et alj[2006l) . The 
observational statistics are dominated by systems in which 
M. > 1O 7 M because kinematic detection of such massive 
objects is easier to achieve. However, if the M-a relation ex- 
tends to lower masses, a possibility supp orted by observations 
of low-luminosity active galactic nuclei dGreene & H0II2OO4I 
iBarth et all2005h . and most correspondingly small spheroids 
harbor MBHs, the nuclear MBHs in the 10 5 - 10 7 M© range, 
of special interest for the present work, may reach a density 
of order 10~ 2 Mpc~ 3 in the local universe jAller & Richstonel 
l200llShankar et all2004l) . 

Our own Milky Way (MW) is the galaxy for which we have 
the strongest observational evidence for the presence of a cen- 
tral MBH. Spectroscopic and astrometric measurements of the 
motion of stars in the vicinity of the radio source Sgr A* in- 
deed indicate that they are orbiting a dark mass concentration 
of some 3-4 x 10 6 M Q whose average density must exceed 
3 x 10 19 M Q pc~ 3 . This high density is incompatible with a 
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stable cluster of an y kno wn less massive astronomical ob- 
jects JC.enzel et aimool: ISchodel et all 120031 iGhez et all 
120051) . and therefore the presence of a central black hole of 
the above mass is well accepted. For a comprehensive review 
of the possible interactions between the Sgr A* MBH and the 
surroundin g stars and their observational consequences, see 
lAlexandedEool) . 

Mass segregation is thought to bring thou sands of stel- 
lar BHs in the innermost pc a round Sgr A* JMorrisl[l993t 
iMiralda-Escude & GoulJ2000l) . This central overpopulation 
may have a variety of consequences. We note that the com- 
pact objects probably dominate the stellar mass density in a 
region (^^OTgcjjn which the "S" stars are confined^Gen- 
zel et al. l2003t ISchodel et alJl2003tlGhez et alJl2005h . From 
infrared photometry and spectroscopy, these stars appear to be 
main sequence (MS) objects with m asses of order 4- 1 6M^ 
and thereforejounger than 100 MvrjGe^arietalj2002LGhez 
et al. T27)03TITElsarihau^r"^^ This apparent youth in a 

environment where normal stellar formation is made impossi- 
ble by the strong tidal forces is an unsolved enigma. 

The compact stars may collide and merge with MS stars 
and giants, hence creating u nusual objects , once suggested 
to be the S-stars themselves JMorrislll993h . or increase the 
rotati on rate of extended stars thr ough multiple tidal interac- 
tions i Alexander & Kumarll2001h . It has also been proposed 
that S-stars are young stars, formed at > 1 pc from the MBH, 
whose orbital eccentricities were increased by some perturba- 
tion such as that of an (unseen) stellar cluster, and which were 
trapped on close orbits arou nd Sgr A* by exchanges w ith less 
massive compact remnants JAlexander & Liviol2004h . 

The presence of the stellar BHs around the MBH can in 
principle be revealed through different kinds of observations. 
If one of the objects acts as a secondary gravitational lens for a 
distant star lensed b y the c entral MB H JChaname et al.l200ll 
lAlexander & Loebl 1200 it I AlexandeJ 120031) . Unfortunately, 
according to these studies, the rate of such double lensing oc- 
curring at a detectable level is very low. 
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If the motions of the S-stars can be tracked with high 
enough a precision, an extended distribution of non-luminous 
matter around the Galactic MBH should signal itself through 
its effect on their orbits. In a slightly non-Keplerian potential, 
the orbits are affected by Newtonian retrogra de precession 
jRubilar & Eckartll200lHWeinberg et alJl2005l) . Present-day 
observations_areJnsufficient to detect this effect (Mouawad 
et al. 120051) . but lWeinberg et*aTT(l2^0^in iave*shown that future 
"extremely large telescopes" (ELTs) with diameters of 30 m 
or more will likely be able to measure the mass and shape of 
a dark density cusp of the form p oc R~ 7 , if it tallies at least 
~ 2OOOM within 10" 2 pc of Sgr A* and has 7 < 2. This 
effect is only sensitive to the overall p(R); it does not distin- 
guish between a population of stellar BHs and another type 
of non-luminous component such as a cold dark matter, al- 
though the latter probably contributes much less than 10 % of 
the density inside the innermost dc of the Galaxv (e.g., Gnedin 
& Primack l2004lBertone & Merrittll2005albl) . However, if a 
concentration of ~ 10M Q BHs is indeed present, ELT obser- 
vations should allow us to witness 2-body relaxation at work 
by the detection of ~ 3 gravitational encounters per year be- 
tween any of ~ 100 monitored S-stars and a stellar BH (Wein- 
berg et al. 120051) . 

Radio pulsars on similar short-period orbits would allow 
the same kind of measurements with a very high accuracy as 
well as precise tests of the theory of general relativity (Cordes 
et al. l2004tlKramer et alJ2004tlPfahl & Loebl2004l). Based on 
the semi-analytical work o f fMiralda-E scude & Gouldl2000l) . 
iChaname & Gouldl J2002I) have suggested that stellar BHs, 
by concentrating around Sgr A* will push out lighter objects, 
possibly creating a central dip in their density profile and have 
pointed out that pulsars would be the ideal probes to detect 
this effect if the sky position of some 50 of them within a few 
arcsec of Sgr A* can be obtained. Unfortunately, because of 
extreme dispersion suffered by radio signals traveling from 
the Galactic center, the detection of pulsars in this region will 
probably require future radio telescopes with hig h sensitivity 
at frequencies > 10 GHz , such as the SKA 4 (ICordes & Laziol 
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H997tlCordes et al.l2004HKramer et al.l200l 

Nevertheless, relatively direct evidence for the presence of 
an abundant population of stellar BHs around Sgr A* may 
not need to await next-generation telescopes. Recently, ob- 
servations with the Chandra X-ray satellite have revealed 7 
transient sources with in 23 pc of projected distance of Sgr A* 
jMuno et al.ll2005bl) : 4 of them have projected distances 
smaller than 1 pc, indicative of an overabundance in this cen- 
tral region by a factor ~ 20 when normalized to the total en- 
closed mass at 1 and 23 pc l lLaunhardt et al.l [20021) . These 
sources are believed to be X-ray binaries, i.e., compact objects 
accreting from a binary companion; however current observa- 
tions do not shed light on whether these are neutron stars or 
black holes with low or high mass companions. However, for 
one case, there is strong evidence for a low-m ass (< 1M^) 
donor and some preference for a BH accret or iBower et alJ 
l200llMuno et all2005atlPorauet et al!2005l) . 

In the present study we undertake a careful numerical 
model exploration of the distribution of compact objects in 
galactic centers. Our goal is to assess the importance and 
detectability of these various effects of mass segregation in 
the context of existing observations. These models are ob- 
tained by explicit integration of the long-term stellar dynam- 
ical evolution of spherical nucleus models with account for 

4 Square Kilometer Array, http : / /www . skatelescope . org 



self-gravity, 2-body relaxation, interactions between stars and 
the MBH, and, in some cases, additional physics such as 
large-angle scatterings, collisions or stellar evolution. The 
models presented here constitute a noticeable improvement 
over the very few si mple estimates of mass se gregatio n avail- 
able i n the literature (lMorrisll 993; Mir alda-Escude & Gouldl 
2000). 

"Extreme-mass-ratio events" (EMREs) in galactic nuclei is 
our other key motivation. EMREs are events where a stellar 
object interacts strongly with an MBH . The best studied case 
so far (first considered bv lHillsl i ll 9751) ) is that of an extended 
star (main sequence or giant) coming so close to the MBH 
that it is partially or totally disrupted by the intense tidal 
forces. The hydrodynamical and stellar dynamical aspects 
of such tidal disruptions have been the object of scores of 
articles (see references at 

http: //obswww.uniq e . ch/~f reitaq/MODEST_WG 4 /TidalD 

for the former aspect and iFrank & Reesl 1197(1 lReesl fr988: 

iMagqrrian & Tremaine|1999nSYer^TjIme^|l999^ Freitag & 

Benz l2002btlWang & Merrittll2004 amongst others, for the 

latter). Although our models also include a simple treatment 

of tidal disruptions and yield rates for these events, we 

are more specifically interested in another class of EMRE, 

namely the coalescence between a compact star and the 

MBH. "Coalescences" as defined here include "plunges" 

when a star suddenly finds itself on a radial, relativistically 

unstable orbit and disappears through the MBH horizon at 

its next periapse passage, and "inspirals" (EMRIs) during 

which the orbit of the stellar object progressively shrinks by 

emission of gravitational waves (GWs) until it plunges. 

EMRIs will be of prime inte rest for LISA ("Laser In ter- 
ferometer Space Antenna", see iDanzman r3 fl996l l2000l and 
http://lisa. jpl.nasa. gov/), the future space borne 
mission to detect GWs with frequencies in the range ~ 10~ 4 - 
0.1 Hz. The waves emitted during the last year of inspiral, as 
the stellar object orbits in the deep gravitational field of the 
MBH, if detected and analyzed successfully, will inform us 
about the geometry of the space time in the immediate vicinity 
of the massive object, thus allowing to probe general relativity 
in the strong field regime, to establish the existence of MBHs 
and measur e with high accuracy their ma sses and spins (iRvanl 
fT9^5llT997llThornell998l:lHughesl200l . 

Predictions of EMRE rates and properties (especially the 
mass of stellar object and the orbital eccentricity when the 
signal starts contributing to the LISA stream) are important 
for the design and of LISA and the development of data anal- 
ysis tools required to extract weak EMRE signals from a 
data stream containi ng noise and a l arge number of other as- 
trophysical sources <Phinne^l27MlPrincdl200llGair et alJ 
120041) . For the GW signal to be in the frequency range of op- 
timal LISA sensitivity during the last year of inspiral (when 
wave amplitude and the interesting strong-field effects are 
the strongest), the MBH mass must be in the range M. ~ 
10 5 - 1O 7 M . In what follows we argue that such MBHs are 
likely to inhabit stellar spheroids in which relaxation time is 
relatively short causing mass segregation close to the MBH. 
This is of great importance for EMREs as the inspiral a stellar 
BH, with a mass ~ 10M Q can be detected in galaxies ~ 10 
times more distant (and therefore ~ 10 3 more numerous) than 
that of a ~ 1 M Q object. 

Determining rates and characteristics of EMREs for LISA 
is beyond the scope of this paper. A few estimates for those 
exist in the literatu re, based on the sam e stellar dynamical 
code as used here (Freitag 2001, 2003) or on other semi- 
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analytical or numerical methods (iHils & BenderlH99l Sig- 
urdsson & Rees}l 997|:|Miralda-Esc ude' & Gouldl2000Hlvanovl 
12001 iHopman & AlexandeJl2005l) . The results from these 
studies are scattered over a disquieting large range, appr ox- 
imately f rom 5 x lO^y r" 1 (IHopman & Alexanderll2005h to 
lO^yr" 1 (lFreitaJ2001l). for EMRIs of stellar BHs in a MW- 
like nucleus (see ISigurdssonll2003l for a brief discussion of 
these various studies and EMRIs in general). This is witness, 
in part, to the lack of realistic agreed-on models for the struc- 
ture of galactic nuclei, causing different authors to adopt dif- 
ferent approximations and values to describe the stellar distri- 
bution around a MBH. With this study, we strive to improve 
this situation. 

Another cause for the disagreement found among EMR1 
studies is the poor understanding of the mechanisms responsi- 
ble for EMRIs. All studies have assumed unperturbed spheri- 
cal galactic nuclei in dynamical equilibrium, in which case 2- 
body relaxation is certainly the main agent for bringing stars 
onto very elongated orbits which may result in EMRIs. At 
th e same time, as already pointed out in the pioneering work 
of lHils & BenderN 19951) and analyzed in detail bv Hopman 
& Alexander encounters with other stars may cause 

an inspiralling star to plunge prematurely, before its orbital 
frequency has entered the LISA band. Although GWs emit- 
ted during the plunge itself will contain some high-frequency 
components, it is unlikely to be detectable by LISA as a re- 
solved source because, typically, tens to hundreds of thou- 
sands of cy cles are required to accumulate enough sign al-to- 
noise r atio llBarack & Cutlerl2004albHGair & Wenl2005l Wen 
& Gair l2005h . Hence, for LISA, the problem is not limited to 
the determination of the rate of coalescences, N coa i; one also 
needs to compute the fraction of those, /emri = A^EMRi/^Vcoai, 
that are "clean" inspirals instead of plunges. However, here, 
we limit ourselves to N coa \, the quantity for which the type of 
simulations we carry out yield robust predictions. We think 
that a real trustworthy estimate of /emri can only be arrived 
at through the use of novel methods, to be developed in the 
future (see §[3}. 

The rest of this paper is organized as follows. In §|2] a quick 
review of the relevant aspects of stellar dynamics in galactic 
nuclei and the previous relevant work is presented. In § [5] 
we describe the numerical method used in our simulations, as 
well as the physics and initial conditions implemented. Our 
main results from our ~ 80 simulations are described in § 0] 
and we conclude in §[5] with a discussion of the astrophysical 
implications of our results and an outlook for future work. 

2. REVIEW OF THE THEORY AND PREVIOUS STUDIES 
2.1. Collisional Dynamics in Galactic Nuclei 

In stellar dynamics, the term "collisional" refers to all situ- 
ations in which the discrete nature of stars, i.e. the fact that a 
stellar system is not composed of a continuous fluid but of in- 
dividual objects, plays a role. In the context of galactic nuclei, 
these effects include 2-body relaxation, direct (hydrodynami- 
cal) collisions between stars and close, dissipative interactions 
between stars and a central MBH. 

In this work, we restrict ourselves to the situation of iso- 
lated, spherical systems in dynamical equilibrium. These as- 
sumptions are made necessary by the numerical method we 
use (the Monte-Carlo code, see § [3} which is still, at the 
present day, the only stellar dynamical scheme able to treat the 
collisional evolution of systems consisting of more than one 
million stars with acceptable realism. They are also adopted 



in almost all other works on the subject because they intro- 
duce a well-defined theoretical framework, allowing in partic- 
ular to make use of methods developed for the study of globu- 
lar clusters. Here we use the term "cluster" for any collisional 
stellar system, including galactic nuclei. Collisional dynam- 
ics, generally with a focus on globular or smaller clusters 
is covered by several textbooks (jBinnev & Tremainelll987t 
ISpitzerlH987tlHeggie & Hutll2003l) : therefore we only recall 
here the few concepts needed to understand the rest of the pa- 
per. More detailed explanations about collisional dynamics in 
the context of M onte-Carlo simulations can be found in our 
previo us papers dFreitag & Benzll20011 l2002bt iFreitag et alJ 
I2OO6I) . 

Barring the effects of mass loss due to stellar evolution, in a 
stationary smooth, spherical potential, stellar orbits would be 
energy- and angular-momentum-conserving rosettes of fixed 
shape and the cluster structure would show no secular evo- 
lution. But the potential is the sum of the contribution of a 
finite number of stars (and a MBH) and is affected by short- 
scale and short-time fluctuations which causes the orbital pa- 
rameters to slowly change. In effect, stars are exchanging 
energy and angular momentum with one another and, to a 
very good approximation this "relaxation" can be idealized 
as due to the sum of a large number of uncorrected 2-body 
encounters leading to small deflection angles. This is the base 
for th e Chandrasekhar theory of relaxation llChandrasekharl 
11960ft on which the Fokker-Planck equation and other approx- 
imate treatment of collisional dynamics are based (Binney & 
Tremaine Tl987HHenonll973h . 

In this picture, one can define a local relaxation time 
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where In A is the Coulomb logarithm, a is the ID velocity 
dispersion, n the number density of stars and (m* ) the aver- 
age stellar mass. The slightly unusual numerical coefficient is 
devised such that a particle of mass (m*) traveling for a time 
St through a field of particles of same mass at a relative veloc- 
ity v re i = V2cr would have its trajectory deflected by an angle 
86, with (59 2 ) = (ir/2) 2 8t/t lU . 

The argument of the Coulomb logarithm is A = £> max /£>o 
where bo = G(m*)/a 2 is the typical impact parameter lead- 
ing to a deflection angle of ir/2 in gravitational encounters 
between stars. In a virialized, self-gravitating system, b max is 
of order the half-mass radius R h and A = with j c ss 0.01 
if stars have a mass spectrum (IFreitag et al.ll2006l and refer- 
ences therein). In the region where the gravitational force is 
dominated by the centr al object, one finds A ss M./(m*) (e.g., 
iBahcall & Wolflfl976l) . In practice, this does not lead to an 
important difference, thanks to the damping effect of the log- 
arithm. For instance, one finds ln(7 c A r *) ~ 15 for N* = 3 x 10 8 
and ln(M./(m*)) ~ 11.5 forM. = 10 5 M o and (m«) = 1M Q . 
Therefore, in most studies, including the present one, a fixed 
value of A (= 7^*) is adopted. C omp arisons with direct 
Af-body integrations, presented in § 14.11 as well as with a 
version of the Monte-Carlo code in which A varies with th e 
distance to the center, from M./(m») to j c N* <Freitagl20 00l 
confirm the validity of this approximation. We set 7 C = 0.01. 

The MBH dominates the gravitational force acting on stars 
within an "influence sphere" with a radius of order GM,a^ 2 
where <jq is the stellar velocity dispersio n at larger distances 
(a more practical definition is given in § 13. 21 for the category 
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of galactic nucleus models considered in our simulations). 
In this central region, the velocity dispersion is Keplerian, 
<j{R) ~ GM m /R. If the stars are distributed according to a 
a power-law density profile, n oc RT 1 , the relaxation time gets 
shorter closer to the MBH when 7 > 1 .5 and longer if 7 < 1 .5. 

In what follows, we call "collision" the event in which two 
stars actually come so close to each other as to touch. Neglect- 
ing deformations due to mutual tidal interactions, a collision 
between stars of radii r\ and corresponds to their centers 
coming within a distance r\ + ri of each other. The cross sec- 
tion for this process is 



5 coii -<n + r 2 ) 
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where m 1 2 are the stellar masses and V K \ their relative velocity 
at a large separation. If field-stars of type "2" have number 
density n-i, and all stars of this type have the same velocity, 
the average time for test-star "1" to collide with one of type 

"2" is 
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To estimate the importance of collisions in the dynamics, we 
assume all stars have the same mass and radius, ra», r* and 
their velocity distribution is Maxwellian with dispe rsion a. 
The collision time is then (iBinnev & Tremainell987l) 
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The numerical relation is valid when the velocities are 
much smaller than the stellar escape velocity, a <C V* = 
{IGmJr*) 1 ! 2 = 617.5 kms" 1 (TO*/M ) 1 < /2 (r*/R©)~ 1/2 so that 
the cross section is dominated by gravitational focusing. This 
ceases to be applicable at distances from the MBH smaller 
than 
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For R < R co \i, er > V*, so the collision time reduces to (note 
the different normalization for n and er) 
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The condition for the collision time to become shorter than 
the relaxation time is also a > V*, for In A sw 10-20. Col- 
lisions at such velocities are unlikely to lead to mergers; a 
fly-by with partial mass loss is the most likelv outcome (Fre- 
itag & Benz l2005l) . Only within R CCI » can collisi ons noticeably 
affect the density profile iFrank & Reeslll976t Sigurdsson & 
Rees 119971) . However, hydrodynamical simulations of colli- 
sions between MS stars show that complete stel lar disruptions 
require a > 5V* and n early head-on geometry (|Benz & Hill; 
fT987l [T99l lLai et al.lH993l iFreitag : & Ben3l2002al l2d05T 
Disruptions are therefore rare and the effect of collisions on 
the stellar distribution is weak, even for^^^£ HJ |^Freitag & 
Benz l2002bl) . 

Gravitational encounters with an impact parameter smaller 
than a few bo lead to deflection angles that are relatively large 
and cannot be accounted for in the standard, "diffusive" the- 
ory of relaxation. Therefore in most approaches, both ana- 
lytical and numerical, these large-angle scatterings have to be 



considered as a separate process. We call them "large-angle 
scatterings" and reserve the word "relaxation" to the effect 
of 2-body encounters with larger impact parameters. On aver- 
age, a star will experience an encounter with impact parameter 
(with /la of order a few) over a timescale 
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The effects of large-angle scatterings on the overall evolu- 
tion of a cl uster are negligible in comparis on with "diffusive" 
relaxation (lHenonlll975t lGoodmanlll983l) . However, unlike 
the latter process, they can produce velocity c hanges strong 
enough to eject stars from an isolated cluster (iHenonll 19601. 
fT969UGTO)drnanl 19831) o r, more importantly, from the "cusp" 
aroundJhecejitral^FiH^jLirij aumg ardt 

et al. I2004al) . Therefore thev may be important for the dy- 
namics of the innermost regions just where mass segregation 
is relevant too. 

A central MBH represents a sink for the stellar system as 
it destroys, captures or -if it forms a very compact binary 
with another object- ejects stars that venture very close to it, 
i.e., within some distance ^i oss . In particular, tidal disrup- 
tion, for a star of radius r* and mass m*, occurs at R\ oss = 
R ti . ~ 1.25 r»(M./m«) 1/3 ( IFreitag & Benzll2002bl and refer- 
ences therein). A quasi parabolic orbit whose Newtonian peri- 
apse distance would be smaller than Ri oss = ^ p iun ae = 8GM.c" 2 
actually plunges directlv through the horizon (Zeldovich & 
Novikov 119991) . Orbits with periapse distance < ^i oss , cor- 
responding to angular momentum (per unit mass) J < /lc — 
y/2GM. Ri OS s, form the "loss cone". For a star with velocity 
v at distance R from the center, the loss cone has an aperture 
angle 9lc = Jlc/(Rv). 

If the star is removed from the cluster in a single close en- 
counter with the MBH, a mature loss cone theory has been 
developed which predicts rates and orbital characteris tics of 
such events JFrank & Reesll976tlB"ahcall & WoUl977t Light- 
man jfeSJianh^jJI^^u^gl^ 

et al. 120041) . The notion of critical radius (R„) is central in 
these cases; it is basically the semi-major axis of an orbit 
for which the relaxation processes cause a change of angu- 
lar momentum per orbital time of order 7lc- Inside R cr loss 
cone orbits are nearly completely depleted ("empty loss cone" 
regime). On the scale of the loss cone, the change of or- 
bital parameters due to relaxation can be treated as a diffu- 
sion process and a direct analogy with the heat equation can 
be used to obtain the average time for a star to be destroyed, 
fdestr.e — ln(^Lc) r rix- At distances larger than R cr , relaxation is 
efficient enough to bring stars into and out of the loss cone 
over an orbital time P or b. The loss cone is therefore full and 
fdestr.f — #LC^orb- The total rate of interactions with the MBH 
is given by V = 4ir J R 2 nt^ stI dR. It peaks around R a for many 
density profiles /?(/?) . 

In cases, such as non-destructive tidal interactions and GW 
emission, in which the star looses energy gradually and is only 
destroyed after a large number of periapse passages, the inter- 
play between relaxation and dissipative processes is not di- 
rectly amenable to the relatively simple loss cone formalism. 
The detaile d analysis of such situations h as only recently been 
pioneered f^jexgn^e^^^ggm^l^SiHopman & Alexan- 
der l2005l) . 

In the sphere of influence of the MBH, orbits of bound stars 
are essentially ellipses precessing on a time scale of order 
(M./M» )0r i,)P rb ^ ^orb where M, is the mass of the central 
object, M^ort, the mass in stars within the apocenter distance 
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of the orbit and P or b the orbital period. On shorter timescales, 
orbits exert torques on each other, thus introducing so-called 
"resonant relaxation" which affects the angular momentum 
on a time scale t Te , ~ (M ,/m*)P mh w In A _1 (M./M* i0r<J )f r i x 
dRauch & Tremainel 19961) . Resonant relaxation is suppressed 
by relativistic precession for very close-by orbits satisfying 
Rperi/Rs < M,/M*, or b with R pa i the periapse distance and 
Rs = 2GM,/c 2 the Schwarzschild radius of the central BH. 
Although resonant relaxation may be much faster than "nor- 
mal" relaxation in the sphere of influence of the MBH, it was 
shown to have only a moderate impact on the rate of tidal 
disruptions in galactic nuclei, because these events are domi- 
nated by st ars with semi-major axi s of order the critical radius 
(see §l3~2l and M,/M *. nrh < 10 jRauch & TremaindH996l 
iRauch & Ingallsll998h . On the other hand, the effects on the 
coalescence of compact objects is likely to be weak due to 
relativistic precession. In any case, the study of this question 
requires a method that can account for non-local gravitational 
interactions between orbits ("2-orbit" effects) and is not un- 
dertaken here. 

2.2. Single-Mass Clusters with a Central Object 

The question of how relaxation will shape the distribution 
of a large number of point-like objects of the same mass or- 
biting a massive object has been addressed in the 70's, shortly 
after the detection of X-ray sources in globular clusters trig- 
ered the hypo thesis that there may be IMBHs at their center 
Peeble1H97l iBahcall & Wo5f!976l LShaoiro & Lightmanl 



gen 

m 



1976; Lightman & Shapiro 1977; Cohn & Kulsrud 1978). The 



approximate solution, first found bv IBahcall & Wor3ll97( 
is this context through a Fokker-Planck-type treatment of the 
stellar dynamics, is the formation of a power-law density, 
n(R) oc RT 1 ^. In this simplified treatment, stars are only de- 
stroyed if they reach a very high binding energy (typically 
£\oss ~ GM m /R t ± for tidal disruptions). As it neglects the 
disruption of stars on (very) elongated orbits (J < /lcX this 
idealized configuration corresponds to an isotropic distribu- 
tion with a zero net diffusive flux of stars in E space and a 
constant outward energy flux 5 . More detailed Fokker-Planck 
treatment accounting for loss-cone effects and other realis- 
tic bounding conditions confirmed the Bahcall-Wolf cusp as a 
very good approximation (Bahcall & Wolf 1977: Lightman & 
Shapiro ll977HCohn & Kulsrudlll978l) . It has since be found 
with other methods also based on the diffusive, l ocal the- 
ory o f relaxation: two types of Monte Carlo codes (|Shapirol 
1985 and references therein; Freitag & Benz 2002b) and a 
gas-dynamical approach dAmaro-Seoane et al.l 120 04;'). Very 
recently, the approximations involved in these computations 
have been vindicated by direct TV-body simulations in which 
the formation of t he R~ 7 / 4 profile over a relaxation time wa s 
indeed witnessed jBaumgardt et all2004atlPreto et al l200l . 

Using a hom ological model for the evolution of a cluster, 
IShapirol i 19771) showed how a central BH can power the ex- 
pansion of the stellar system, by destroying stars that have 
diffused deep into the cusp. A central BH is therefore able 
to drive gravothermal expansion in a way similar to harden- 
ing binaries, but without leading to core oscillations CHeggie 
& Hut 120031) . The central BH acts as a heat source for the 

5 Treating the cluster as a conducting gas, the Bahcall-Wolf solution can 
be found by imposing dF /dR = with F = — 4ttR 2 ■ K(da 2 / dR) the rate of 
"thermal" energy conducted across a sphere of radius R. re is the thermal 
conductivity, re = pA 2 /r with p the mass density, A R the effective mean 
free path and r ~ f r ix the timescale for energy exchange. 



whole cluster only if, on average, destroyed stars have neg- 
ative orbital energies rel ative to the BH, a con dition roughly 
equivalent to R cl < R^ dDuncan & Shapiroll98l . 

2.3. Multi-Mass Clusters with a Central Object 

Surprisingly, the effects of relaxation in a multi-mass clus- 
ter containing a central massive object have been little stud- 
ied. To our knowledge, the only in-depth theoretical study of 
mass se gregation in the Keplerian potential of a MBH is the 
work of B ahcall & "V^ll ( fT977l) . The long-term evolution of a 
few models of MBH-hosting galactic nuclei with a mass spec- 
tru m was followed num erically using a Fokker-Planck code 
bv lMurphveTall i ll 9911) and with the same Monte-Carlo code 
as the present study by iFreitag & Benzl J2002bh and iFreitad 
However, in those studies, rich physics was included 
which complicates the interpretation of the results (collisions, 
stellar evolution,. . . ) and their authors did not present detailed 
re sults concernin g mass segregation. Also, with the exception 
of lFreitaglll200l . the initial conditions used were not tailored 
to represent anv specific galactic nucleus. Recently, (Baum- 
gardt et al. I2004hl) carried out direct TV-body simulations of 
multi-mass clusters with some 1.6 x 10 4 to 1.3 x 10 s stars 
hosting a central IMBH and discussed how stars of different 
masses distribute themselves around the central object. This 
study offers the only direct characterization of mass segrega- 
tion around a massive object. One should be cautious, how- 
ever when trying to apply these TV-body results to larger sys- 
tems such as galactic nuclei because small-TV effects (large- 
angle scatterings, binary interac tions, IMBH wanderin g,. . .) 
may play a significant role there ILin & TremaineH"980h . 

A first step towards the understanding of mass segregation 
in galactic nuclei is to consider the simpler problem of the 
evolution of one or a few massive "tracers" in a non-evolving 
stellar background. We undertake this step here for illustrative 
purposes. This is a useful idealization for the early dynamical 
evolution of the population of stellar BHs. Those are very 
rare objects so, until they have concentrated in the innermost 
regions, they will mostly interact with other stars and not with 
one another. We assume all stellar BHs have mass m^n and 
all other stars have mass m with q = otbh/'« ^ L 1 = 30 is a 
realistic value. The effects of 2-body relaxation on the orbit 
of a massive particle ("test particle") in a field of much lighter 
field particles is e mbodied in the classical dy namical friction 
(DF) formula (see lBinnev & Tremainel 19871 § 7.1) 

., _ 47rln AG 2 nm(m + mBH) 



a DF -- 



'DF' 



-k(X)v (8) 



with k(X) = eTf(X)-2TT-^ 2 Xe- xl and X = v/(V2a). In this 
formula, adf is the force per unit mass on the test particle 
due to DF, v its velocity, n the density of field particles, and 
a the (ID) dispersion of their velocities, assumed to have a 
Maxwellian distribution. f DF is of order the local relaxation 
time divided by 1 + q so the massive particles should already 
experience significant mass segregation after a small fraction 
of the relaxation time. 

For an object on a circular orbit of radius R, v = v c = 
\J GM enc i(R) I R where M enc \(R) is the total mass within R, and 
a differential equation for the evolution of R is easily derived 
from «df^ = d(Rv) /dt, 

dR_ ( 2nGRhi(R) m V 1 R 
dt \ v c (R) 2 t DF (R) W 

Although it can yield a qualitative understanding and a 
first approximation to the development of mass segregation, 
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a treatment based on the use of Eq.|9] falls short of >hysical 
realism. First, relaxation reduces to dynamical friction only 
in the limit of very large mass ratio. In general, the c irection 
of v (and not only its modulus) is also affected by 2-1 ody en- 
counters, causing the eccentricity of a circular orbit to drift 
away from zero. Second, if massive objects are ni merous 
enough, they will eventually come to dominate the central re- 
gion. There, they will push the lighter objects away by heating 
them and start interacting with each other in a way mpre sim- 
ilar to the single-mass situation. The dynamical fricf on pic 
ture does not provide a way to determine the quasi-stationary 
distribution the particles of different masses will adof t on the 
long term. 

In another seminal paper, iBahcall & WoiH ( fT977i) | 



studied 
the po- 
steady- 



the possibility for a multimass system dominated by 
tential of a central MBH to settle into a relaxational 
state configuration (a cusp), provided stars lost to interactions 
with the MBH are replaced by stars coming from rr ore dis- 
tant regions. By solving the coupled Boltzmann equa ions for 
stars of various masses, they found that the stars o ? differ- 
ent masses, m,-, should approximately follow one-part icle dis- 
tribution functions that are power law of the binding energy, 
fi(x, v) = F/(E) oc E Pt , with indices scaling like 



El 

111. 



El 

m; 



(10) 



7=1 + 
ts, who 
single- 



These correspond to density profiles «, oc R~ J with 
p. They found p ~ 0.30 for the most massive obje< 
dominate the central density, close to the value for 
mass distribution, p ~ 0.25. For much lighter objects in th e 
innermost regions, 7 ~ 1 .5 is expected (see also Merrijtl2004l) . 

It is interesting to note that the massive stars concentrate 
to the center because they lose energy to lighter ones dur- 
ing 2-body encounters. This tendency would yield statisti- 
cal equipartition of kinetic energy if it wasn't for the overall 
gravitational potential in which the heavy objects sink, thus 
increasing their velocities. In a cluster without a cenrnal black 
hole, equipartition can only be reached at the center end only 
if the massive particles are in small number or have a r lass not 
much exceeding that of the lighter ones, so that the) cannot 
form a self-g ravitating system, with negative heat c apacity, 
on their own (lSpitzerll969t | Vishniadl978tlln"agaki 8a Saslawl 



1985 



IWatters et alJl2000UGurkan et al.ll2004HKhal 



si et al 



2005). For all realistic mass spectrum, mass segrega ion will 
trigger the core collapse of the sub-system of massive bodies, 
a process known as "Spitzer instability" 

Clearly, in a fixed Keplerian potential, massive stars can 
never reach equipartition with lighter ones; as they concen- 
trate to the center, their velocity dispersion must ncrease 
and the thermal imbalance with the lighter objects is main- 
tained. An accelerated, catastrophic collapse of the popula- 
tion of massive objects is prevented, however, by the heating 
effect of the central MBH which eventually compen? ates for 
the energy lost to the light stars. Hence, a cusp of massive ob- 
jects is expected to form and maintain itself in therrr al quasi 
equilibrium while it drives the expansion of the distribution of 
lighter objects. 

Published simulations of multi-mass clusters with $ central 



I)MBH are few and far apart. The work of iMurp] iv et al.l 



1991) stands out as a pioneering effort to follow th? evolu- 
tion of galactic nuclei taking into account relaxatior , stellar 
evolution and collisions. These authors have published lim- 
ited data from one run without stellar evolution or collisions. 



& Wolfl JI9771) relative to the cusp exponents for stars of dif- 
ferent n lasses (Eg. HOt . From their Figure 9, however, it seems 
that the region for which this applies encompasses of order 
10 4 M G only, at a time when, judging from their case 4C, the 

MBH has certainly g rown past 10 6 Mp . 

To PL r knowledge. lBaumgardt et al.l J2004bh have presented 
the onl ' direct A/-body simulations of a multi-mass system 
with a central massive object. Although they observe that 
the mo; t massive objects form a power-law cusp of exponent 
compatible with 7 = 1.75, the central profiles of the lighter 
species are found to be much shallower than predicted by 
Ea.lTol with 7 ~ J5 + m*/(lAMn). Ho wever, in the light 
of our comment on iMurohv et alt ill 99 11) and of our simu- 
lationsjtlriscannot be interpreted as a rebuttal of Bahcall & 
Wolf ( 1)9771) but more likely is an indication that the appro- 
priate regime is only reached deep in the influence region, a 
region not probed by A/-body simulations with N p < 131 000. 

SIMULATIONS: METHOD AND INITIAL CONDITIONS 



This 



The Monte-Carlo Code for Nucleus Dynamics 



They report a good agreement with the prediction of lBahcall 



oil 
f B 



work is based on simulations of the long-term stel- 
lar dynamical evolution of galactic nuclei performed with 
ME(S5 Y)**2. This code is based on the Monte-C arlo al- 
gorithm first described and implem ented by Henon l lHenonl 
11971 blit iHenorj 1 973|: |He nonl 1 9751) . It has been described in 
detail b ^ iFreitag & BenzH2001tl2002bl) . Here, we succinctly 
remind the basics of the method and the included physics. 

The Monte-Carlo method is based on the assumptions of 
spherical symmetry and dynamical equilibrium. The cluster 
is repre rented by a number (typically N p = 10 5 - 10 7 ) of par- 
ticles, e ach of which is a spherical shell. These shells consti- 
tute a sampling of 1 -particle the distribution function in the 
phase a nd stellar-parameters spaces. In other words, a shell 
corresponds to stars with a given orbital energy E and an- 
gular rr omentum (in modulus) J and given stellar properties 
(mass, ige, etc.). At any time a shell also has a given radius 
R. Eacl shell represents the same number of stars, N^/Np > 1 
(for small systems, one may set N p = N* and N p > N* is for- 
mally possible). 

Orbital motion is not followed as dynamical equilibrium is 
assumed (the system is phase-mixed); instead, the position of 
a particle on its orbit, i.e. its radius R, is selected with proba- 
bility reflecting the time spent at each R on the orbit specified 
by E and J in the potential of the other shells and the central 
object. 

Grav tational relaxation is treated in the Chandrasekhar pic- 
ture, similarly to what i s done to derive the orbit -averaged 
Fokker- Planck equation dBinnev & Tremainell987l) . It is as- 
sumed to reduce to the effect of a large number of uncor- 
rected, small-angle 2-body scatterings dominated by impact 
p aram eters bo <C b <C £> max (the value of b max is discussed in 
§ 12. 1> . Consequently, relaxation is implemented as a series 
;ity perturbations between neighboring particles. In 
Y)**2, time steps 5t are a function of the radius R 
set to be smaller or equal to a fraction fs t of the local 
f r ix. F01 the present work, we set f$, = 0.04 and checked that 
fst = 0. )1 does not lead to significantly different results. 

Stelkr collisions can also be treated by computing the col- 
lision time for a pair (Eq. [3} and comparing 8t /t co \\ to a uni- 
form [Q, 1[ variate. For the simulations of the present work 
where collisions were included, interpolation from a large 
databas e of Smoothed Particl e Hydrodynamics (SPH) simu- 
lations iFreitag & Benzl2005l) was used to determine the out- 
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c ome, as described inlFreitag et alJ (12006ft . The simulations 
of lFreitag & Beiizl J2005h" specifically probe the high- velocity 
regime found in the vicinity of MBHs. 

An accurate treatment of the loss-cone process is not possi- 
ble in the framework of the present version of ME(SSY)**2 
because it would require to endow particles in on near the 
loss cone (or on orbits eccentric enough to possibly lead 
to EMRIs) with time steps shorter than the timescale taken 
by relaxation to modify significantly the pericenter distance 
fr.p — (1 — e)f r ix <C frix- St(R) must be an increasing function 
so that setting a short time step for some particle would, in 
practice, reduce the time steps of all particles with positions 
lying inside its apocenter. This difficulty is circumvented by 
an approximate treatment of the relaxation-induced random- 
walk of t he direction of a particl e's velocity vector during a 
time step JFreitag & Benzl2002bl) . 

A novelty introduced in a few runs presented here is the 
treatment of large-angle scatterings. They are treated in a way 
similar to collisions but with a cross-section 

, (1 , 2) _ G(mi+m 2 ) 



it u( l - 2 )\2 
= TT(/lA&o ) 



with = 



Vz 



(11) 



When a large-angle scattering is deemed to occur, the im- 
pact parameter b is selected at random between and /la&o 
with probability density dP/db oc b. The outcome, in the 
center-of-mass frame of the pair, is a deflection of the veloc- 
ity vectors by an angle 2arctan(Z?o ,T) jb). When large-angle 
scatterings are included, the Coulomb logarithm is reduced to 
ln(7 c A^ //la) to account for the fact that gravitational encoun- 
ters with b < /la^o are now treated separately. 

3.2. Initial Nucleus Models 

As is customary in c luster simulations, we use the 
'W— b ody" system of units ( !Henonll971atlHeggie & Mathieul 
119861) . Unlike the situations for which this system was first in- 
troduced, we deal here with stellar systems that are not strictly 
self-gravitating; instead their central regions are dominated by 
the potential of a massive, fixed object. Hence we define the 
unit system such that the constant of gravity is G = 1, the total 
stellar mass is initially M c \(0) = 1, and the total initial stel- 
lar gravitational energy (not accounting for the contribution 
of the MBH to the potential) is -1 /2. We denote by R NB the 
A/-body length unit. 

As a time unit, we use the "Fokker-Planck time" Tpp which 
is connected to the A^-body time unit Tnb through Tpp = 
(A/*(0)/ln A)7nb were Af*(0) is the initial number of stars. We 
prefer to use Tpp rather than Tnb because the former is a relax- 
ation time while the latter is a dynamical time. We consider 
systems in dynamical equilibrium whose evolution is secular, 
in most cases driven by 2-body relaxation. For a large vari- 
ety of cluster str uctures ? Tpp 10f r h where f r h is the half-mass 
relaxation time JSpitzerll987T) . 



frh = 



0.138AT* / Rl 



1/2 



(12) 



In A VGMci, 
with Rh the radius enclosing half of the stellar mass. 

There are only few published models for (spherical) clus- 
ters in dynamical equilibrium and containing a massive cen- 
tral object. The best described and most conve nient ones are 
the "eta-models" introduced bv lDehne n|(|1993l) and e xtended 
to systems with a central object bv iTremaine et alJ d!994l) . 
The density profile is 

»7-3 / D \ — TJ— 1 



p(R) = 



1 + - 



R 



(13) 
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Fig. 1 . — Influence radius R m f[ as a function of the parameters fj, in eta- 
models. Each curve is for a value of r\. The dots indicate the value of the 
break radius R^. This diagram allows to find the maximum fi value for which 
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FIG. 2. — Enclosed mass as function of radius for some of our models. The 
points in gray are observational constrains from the kinematics of stars and 
gas at the center of the MW. The point at the smallest distance corresponds 
to the simultaneous fit of the orbi ts of S-sta rs by Ghez et al. 1 2005). Other 
data points have been compiled by Schodel et al. 1 2003). The thin lines show 
the stellar contribution and thick lines include the central MBH with M, = 
3.5 X 1O 6 M . The solid line is our reference model. The short-dashed line 
represents a model with a total stellar mass two time smaller but same stellar 
density at small distances (R <C Snb)- 



The exponent 77 can take any value between and 5/2. At 
small radii, p oc R' 1 with 7 = 3 - 77 while at large distances, 
density falls off like R~ A . The break radius can easily be ex- 
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FIG. 3. — Time scales in our reference model of the MW nucleus (M c i = 
7 X 10 7 Mg, r) = 1.5, i?NB = 28 pc). We plot the orbital time f or t, (assuming 
a circular orbit), the relaxation time f r i x , the time scale for large-angle de- 
flections f] a , collision times f co n and the timescale for diffusion by 2-body 
relaxation over the loss cone for tidal disruption, fLc = $Lc f rlx • ^ or c °Hi s i° ns 
we indicate the average time for a MS star to collide with another MS star or 
with a red giant (RG) and the average time for a RG to collide with another 
RG. We assume M MS = M RG = 1 M , R MS = 1 R , Rrg = 50M Q , and 5 % 
of RGs. The radius where f or b = fLC is the critical radius. 



pressed in terms of other important length scales, 



R h = (277- 1) R 



NB 



:(2 



(14) 



The fraction of the stellar mass enclosed by R\, is 2 i. The 
central MBH defines a second dimensionless parameter /i = 
M./M ch 

At short distances from it, the MBH dominates the dynam- 
ics and therefore a 2 (R) w cr^ BH (/?) = {4-rjT l GM./R. We 
define the influence radius R m f[ implicitly through u 2 (/?i n fi) = 
2cr^ IBH (7?i n fl). Figure^shows how R m ^ depends on fx for var- 
ious values of 77. In the present study, we use eta-models as 
a way to carry out simulations with a power-law density cusp 
of controlled exponent 7 as initial conditions. We view the 
steeper density decrease at large radii, R > R\,, as a cut-off to 
avoid wasting computer memory and CPU time by putting a 
large number of particles at distances that should not be influ- 
enced by the presence of the MBH through relaxation effects. 
In other words, the value of R\, should be irrelevant as long 
as it is large enough to encompass the region within which 
the collisional physics takes place. It is therefore important 
to have R\, > Ri n a and, from Figure [2 we see that this will be 
the case for 77 = 1-2.25 provided that /.t < 0.05. For 77 < 1.5, 
yU < 0.1 should be sufficient. 

Another impo rtant radius is the critical radius for tidal dis- 
ruptions. R,- r ,j jFrank & Reesl 119761 iLightman & Shapiro! 
T977t iMagorrian & Tremaindll999t fSver & Ulmerlll999t 
Amaro-Seoane et afll2004l) . It is defined as the position in 
the cluster where the diffusion angle caused by relaxation per 
orbital time equals the opening angle of the loss cone, for a 



(1) 



tibt 



- U LC- 



(15) 



A local, typical value of #lc can be obtained by computing it 
for a star whose velocity would be equal to the (3D) vel ocity 
dispersion, i.e., solving Eq. 23 o flFreitag & BenzH2002bl) with 
v 2 = 3a 2 (R). 

The rate of tidal disruptions is dominated by the contribu- 
tion of stars with apocenter distances of order of the mini- 
mum between R m f[ and R aX d- For all models considered in 
this study, R aX d < Rinft < Rb (see Figure [3ji so the loss-cone 
effects should be little affected by the existence of a steeper 
density decrease beyond R\,. 

For the present study we construct most models so that 
they best approximate the conditions in the MW nucleus. In 
Figure 13 we plot the enclosed mass as a function of radius 
for some o f our initial models and compare wi th observa- 
tional data JSchodel et all2003tlGhez et alJl2005h . Our refer- 
ence cluster model is described by M. = 3.5 x 10 6 M Q , M c \ = 
1 x 10 7 M Q (hence \x = 0.05), 77 = 1.5 and R NB = 28pc. This 
model has a central density cusp with 7 = 1.5, a va lue consis- 
tent with the stellar coun ts at the galactic center ( lAlexandeJ 
ll999UGenzel et alJl2003T) . However, a detailed modeling of 
the Galactic center is not our goal. This would in particular re- 
quire ad hoc assumptions regarding the history and locations 
of sta r formation to acc ount for a population with a variety of 
ages dFiger et alJ2004t) . This variety could be the result of the 
intermittent formation, at a few pc from the center, of small 
clusters that th en spiral in and deposit their stars in the nu- 
cleu s (see. e.g..lMaillard et alJ2Q04HPaumard et al.ll2004t Lu 
et al. 120051 for observations andjKim & Morrisll2003t Porte- 
gies Zwart et al. l2003tlKim et alJ2004HGurkan & Rasiol2005l 
for simulations). 

Separate from the MW-like models, we explore the effects 
of mass segregation in idealized galactic nucleus models with 
a variety of structural parameters. We consider nuclei with 
M. in the range 10 5 - 10 7 M Q . To decrease the dimensionality 
of the parameter space, we assume the M-a relation (Merritt 
& Ferrarese l2001tlTremaine et al.l2002tlBarth et alJl2005l) to 
hold perfectly, 



M. ~M 



100 



lOOkms" 



(16) 



Tremaine et al. <2002l) find /3 = 4.0 ± 0.3 and M m = 
8.3*33 x 10 6 M©- With a velocity dispersion of a ~ 
lOOkms -1 , the MW harbors an under-massive MBH. We 
note, however, that the velocity dispersion of the MW cen- 
tral region, as defined for use in Eq.^J is do minated by stars 
locate d at a few hundreds pc from the center llTremaine et alJ 
I2002L and references therein), a region we do not attempt to 
model. The MW nucleus is the only one whose structure is 
relatively well constrained by observations at the scales of in- 
terest here (R < lOpc). Hence, for simplicity we adopt the 
MW nucleus as typical. A model for a nucleus with an MBH 
of mass M. is obtained from a MW model of same 77 and /j, by 
simple length (and mass) rescaling. As R oc Ma~ 2 and using 
(3 = 4 we obtain: 



^nb - ^nb| mw 



M. 



3.5 x 1O 6 M 



1/2 



(17) 



where /?nbImw I s A^-body radius of the MW model. 
Neglecting the dependence of the Coulomb logarithm on 
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N*, the relaxation time of the model scales like f r j x = 
tru\ MW (M./3.5 x 10 6 M o ) 5 / 4 . It follows that t dx (R m ) ex- 
ceeds ~ lOGyr for M. > 10 7 M Q and we expect only minor 
relaxation effects in such massive nuclei. The lowest MBH 
mass we consider, 10 5 M Q , corresponds to that of the small- 
est M BH detected with some c onfidence in a galactic center 
so far dBarth et al.ll2004ll2005l) . Even smaller systems, such 
as the nuclei of dwarf galaxies or globular clusters may host 
intermediate mass black holes (IMBHs) with M, < 10 5 M Q . 
We do not address the evolution of low-mass objects here be- 
cause their central dynamics may be significantly influenced 
by small-Af* effects not included in ME(SSY)**2. Recently, 
the very significant increase of comput ational power offere d 
by special-purpose GRAPE hardware (Makino e t al J 12003b . 
combined with a variety of mathe matical and num erical tech- 
niques to speed up computations jAarsefhll2003l) have made 
it possible to follow the relaxation evolution of clusters with 
a central massive obj ect and up to 2.5 x 10 s stars by "d irect" 
Af-bodv integration s jBaumgardt et al.ll2004albL l2005t Preto 
et al. 12004) . However, because of the steep dependence of 
the CPU time on the number of particles imposed by direct 
force computations in the A^-body algorithm (approximately 
?cpu oc Afp per relaxation time), systems containing ~ 10 6 stars 
or more can only be studied with more approximate methods 
such as MC codes. 

The range M. = 1 5 - 1 7 M also corresponds to the MBH 
around which an EMR1 has the best chance to be detected by 
LISA. The orbital period of a test particle on the innermost 
stable circular orbit around a non-rotating BH of mass M. is 

c 3 / M V 1 

^^2^^r 2 - 2xl0 ~ 3Hz {wk) ■ (18) 

Consequently, inspirals into MBH more massive than ~ 
1O 7 M produce signals with frequency too low for LISA 
to detect, while the final inspiral into an MBH with M. < 
10 5 M Q occurs at periods higher than the time taken by light 
to travel along LISA arms, which strongly reduces sensitivity 
at those frequencies tLarson et alJl2000» . In principle EM- 
RIs into such lower-mass MBH could be caught at an earlier 
phase in their orbital evolution but the emitted waves have 
much lowe r amplitude then, thus severely limiting the detec- 
tion range l|Wilj200j| i F^urthermore, the anal^sis^crfHopman 
& Alexander (120051) indicates that most stellar objects closely 
bound to an IMBH will be scattered on to a direct plunge orbit 
before they enter the LISA band. These authors predict that 
successful (i.e. gradual) LISA inspirals around IMBH have to 
start at very high eccentricities and small semi-major axis and 
should last only ~ 1 year before coalescence. 

In this study, we concern ourselves with the idealized situa- 
tion of an isolated, gas-free galactic nucleus. In particular, we 
do not consider the effects of interactions with other galaxies, 
such as mergers with other nuclei or gas inflow. Similarly, 
we neglect the possibility of smaller stellar clusters spiraling 
down to the galactic center or non-spherical mass distribu- 
tions. Finally, we assume that all stars have formed in a sin- 
gle burst with no further star formation. For the models in 
which stellar evolution or collisions are included, the gas lost 
by stars is considered instantaneously lost from the system, 
with, in some cases, a fraction being accreted by the central 
MBH. 

Some of these simplifications, most noticeably that of 
spherical symmetry and absence of gas, are required by the 
numerical methods used. Others are made in order to reduce 
the complexity of the problem and the dimensionality of the 



parameter-space, hence allowing a better understanding of the 
systems under study. 

Most of our simplifying assumptions favor mass segrega- 
tion of stellar BHs, our primary object of study. For instance, 
it seems likely that a merger between nuclei induces vio- 
lent relaxation, thus erasing -at least partially- any previous 
mass segregation. If both nuclei contain a MBH, the binary 
MBH will eject stars from the central regions and strongly 
decr ease the density there, thus lengthening relaxation tim e 
(e.g.. lMilosavlievic & Merritl200ltlMakino > & Funatol2004l) . 
Also, if stars are formed over an extended period of time in- 
stead of all being born at some "initial" time, stellar BHs will, 
on average, have less time to experience mass segregation. 

Cosmological simulations indicate that most normal galax- 
ies have not suffered a major merger for several Gyrs. In 
particular, some 5-7 Gy r are probably required for a disk to 
(re)fo rm after a merger (iGovernato et al.ll 19941: lAbadi et alJ 
120031) . Therefore our simulations can be considered to cover 
the evolution of a galactic nucleus since it experienced its last 
major merger. We will focus our analysis on the structure of 
the nucleus after 5 and 10 Gyr of simulated evolution; 5 Gyr is 
a reasonable value for the period of time during which a nu- 
cleus in the present-day universe may have evolved without 
strong interactions; 10 Gyr is an upper limit that enables us to 
see what the maximum effects of relaxation are likely to be. 
Mergers probably lead to important gas inflow into the cen- 
tral regions, triggering stellar format ion and accretion onto the 
MBH, in a complex interplay (e.g.. lSpringel et a!1l2005l) . In 
such episodes the MBH may grow substantially on time scales 
shorter than the relaxation time, but still significantly longer 
than stellar orbital periods. The stellar nucleus then contracts 
adiabaticall y in response to the deepening of t he MBH poten - 
tial (lYoungll980HOuinlan et al.ll995UFreitag & Benzl2002bl) . 
To investigate the impact of such episodes on the structure of 
the nucleus several Gyrs later, and contrast it with our stan- 
dard models where the mass of the MBH increases only little 
during the course of the simulation (by tidally disrupting and 
capturing stars), we computed a few models in which a central 
BH of small mass (/i(f = 0) = 10~ 5 ) grows rapidly by accreting 
some fraction of the gas released due to stellar evolution. 

3.3. Stellar Population and Evolution 

Except for a few test-case models presented in § 14.11 we 
use the "Kroupa" initial mass function (IMF ) for all our mod- 
els JKroupa et al.lll993t iKroupal 1200 laid) . It is a broken 
power-law, dN*/dm* oc m~ a , with a = 0.3, 1.3 and 2.3 in the 
ranges m*/M s G [0.01,0.08[, [0.08, 0.5[ and [0.5,120], re- 
spectively. We generally consider the range 0.2- 120 M for 
stellar masses on the MS. 

In most simulations, we do not include stellar evolution but 
start with a stellar population in which all stars already have 
an age of 10 Gyr. This is of course not a physically consis- 
tent treatment but we choose it for the sake of simplicity. For 
comparison purposes, in a few simulations, stellar evolution is 
included and those simulations are started with zero-age MS 
stars. The main impact of stellar evolution is to induce signif- 
icant mass loss in the first ~ 10 s years. As we will see, the 
nucleus experiences strong expansion if this gas is expelled 
from it. To produce such a model for a nucleus with specific 
current properties (as those of the MW), we have to find, by 
trial-and-error, initial cluster structural parameters leading, af- 
ter 5 - 10 Gyr, to a nucleus model fitting the observations (in 
our case, the enclosed mass as function of radius). 

We use a simple stellar evolution prescription according to 
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which stars keep a fixed mass and radius while on the MS 
and instantaneously turn into compact remnants (CRs) at the 
end of their MS lifetime, tusim*)- Data for fmsQw*) were pro- 
vided by K. Belczynski l lHurlev et alJl2000uBelczvnski et alJ 
As for the relation between the stellar mass on the MS 
and the nature and mass of the CR, we consider three models, 
presented in Tabled In the first one, dubbed "fiducial" (F), 
we assume all white dwarfs (WDs), neutron stars (NSs) and 
stellar BHs have a mass of 0.6, 1.4 and 10M Q respectively. 
Th e two other model s make use of the prescriptions developed 
bv lBelczvnski et all (120021) . assuming either solar (Z = 0.02, 
model "BS") or metal -poor (Z = 10 -4 , model BP) chemical 
composition. These prescriptions represent our current under- 
standing (although incomplete) of massive-star core collapse 
and possible fallback onto the nascent compact remnant. The 
quantitative aspects a re consistent w ith the hydrodynamic cal- 
culations presented in lFrveJ i 19991) and the resulting relations 
between MS and CR mas^es_are ii showriJri i J'igure 1 of_Bel- 
czynski et al. i2002ft . When stellar evolution is included, we 
impose that the time step is smaller than a factor /se times 
the MS life time t M s for all stars still on the MS. We have 
set /se = 0.1 after have checked that results are essentially the 
same as with /se = 0.025. 

To explore the effect of supernova kicks in some simula- 
tions with stellar evolution, we give NSs and BHs a veloc- 
ity kick at birth. Although the mechanism responsible for 
such "natal kicks" is still not understood, they are required 
to explain the high spacial velocities of observed field pulsars 
l lHobbs et al.l2005L and references therein) as well as other ob- 
served characteristics of neutron star binaries (e.g., Willems 
et al. l2004tlThorsett et al.l2005L and references therein - ). 

There are also observations and interpretatio n analyses sug- 
gesting that some BHs receive a kick at birth jMjrab el et alJ 
I2TO 12001 IGualandris et aTll2()05aLlWillems et aLKOOlH t 
is generally accepted that a supernova explosion is required 
to provide the natal kick. Consequently, it is likely that only 
BHs formed through the fallback mechanism, with a progen- 
itor less massive thari^2£gW^2M ai jgc^iy^^icks ii (Fryer & 
Kalogera T200lUHeger et al.ll2003l) . In the MC simulations 
with natal kicks, w e base our prescription on the results of 
iHobbs et al.H2005l) . The kick velocity is picked from a single 
Maxwellian distribution with a one-dimensional dispersion of 
cjnk = 265(1. 4M /m)kms _1 where m is the mass of the NS 
or BH. BHs resulting from the evolution of a MS star more 
massive than 42 M are not given any kick. 

4. RESULTS OF SIMULATIONS 

Our simulations fall into two categories. First are a few 
cases with a single-mass or a two-component stellar popula- 
tion. They are used to test the MC algorithm by comparison 
with analytical or A^-body results. The second category con- 
sists of more than 80 galactic nucleus models with more re- 
alistic choices of parameters and stellar populations. In what 
follows we describe the results of some representative runs 
and explain how the important outcomes are affected by the 
initial conditions and physics. 

4.1. Test Models 

Since ME(SSY)**2 was or iginally developed and tested 
JFreitag & Benzll200H l2002bh . the code has gone through 
many small revisions. Furthermore, at that time, only few 
direct A^-body simulations had been published with high 
enough resolution to yield test cases to which the results of 
the more approximate MC code could be usefully compared. 
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FIG. 4. — Evolution of the density profile for a single-mass cluster with 
r\ = 2.25 and fj, = M,/M c [ = 0.05 simulated with 4 million particles. Note the 
establishment of a p oc R -175 cusp and the expansion of the cluster, driven by 
diffusion of stars towards the MBH. 




Time [FP units] 

FIG. 5. — Comparative evolution of single-mass models with different ini- 
tial central density profiles. We plot the Lagrange radii, i.e. the radius of 
spheres enclosing the indicated fraction of the total stellar mass. Models with 
77 = 1.5 (solid lines) and r\ = 2.25 (dash-dot) are compared. Both runs have 
[i = 0.05 and N p = 4 X 10 6 . At late times, the two cases have converged to the 
same structure and evolution. 



The advent and spectacular increase of computing speed of 
GRAPE boards now permits more comparisons, although re- 
strictions in the applicability of comparisons still exist. We 
have recently carried out new tests for the core-collapse evolu- 
tion of cluster s with a variety of s tellar mass functions but no 
central object (tFreitag et al.l2006T) . Here, we investigate mod- 
els with a central MBH. We compare MC results with simple 
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FIG. 6. — C omparison of the Lag range radii evolution between /V-body 
results from Baumeardt et al. 1 2004a) and our ME(SSY)**2 simulations of 
single-mass cluster models. Fr om top to bottom, the pan els correspond to 
models 1, 2 and 4 in Table 1 of Baumeardt et al. 1 2004a). The MC results 
are plotted with solid lines, the A'-body results with dashed lines. We present 
MC results obtained with N f =N*= 80000 for case 1 and 4 (top and bottom) 
and N p = 5N« = 320000 for case 2. A run with N p = N* = 80000 gives 
very similar (but noisier) results. The N— body time unit is converted into 
FP unit assuming 7 C =0.11. The curves have been smoothed with a running 
averaging procedure using a Gaussian Kernel. 



semi analytical predictions as well as published and original 
Af-body simulations, presented here for the first time. 

The development of a p oc R~ l 75 density cusp in a single- 
mass cluster hosting a central MBH has bee n a well-accepted 
theore tical prediction for nearly 30 years dBahcall & Wolfl 
11976ft . but h as only recently been verified by direct jV-b ody 
simulations jPreto et al1l2004t iBaumgardt et alJl2Q04ah . In 
Figure |4] we show how such a profile forms in one of our 
single-mass MC simulations of a cluster model with r\ = 2.25, 
fi = 0.05 and jV p =4x 10 6 . It is evident that at late times, 
the evolution is an approximately self-similar expansion of 
the cluster, driven by destruction of stars by the MBH (whose 
mass was kept constant in this simulation). Models with dif- 
ferent initial rj values converge to the same structure and evo- 
lution after t ps (0.05-0.1)^, as illustrated in Figure[5] To 
measure the speed at which the central regions evolve, the re- 
laxation time at the influence radius, U\ x {R m a) (using Eq.|2j(, is 
a more relevant timescale than 7pp; we find f r i x (^inn) = 2.2 x 
10" 3 r FP for (j, = 0.05, ri = 1.5 and f r i x (7? infl ) = 4.9 x 10" 3 r FP for 
77 = 2.25. Hence, the full development of a Bahcall-Wolf cusp 
requires of order lOfi-ixC Rjnfi) in a single-mass clus ter. 

With an A^-body code. Baumgardt et alJ l l2004al) have com- 
puted the evolution of single-mass Wn = 10 King models (Bin- 
ney & Tremaine n*987tlHeggie & Hutl 120031) with a central 
BH of fi in the range 0.0026-0.1 (the stellar velocities were 
modified to ensure approximate dynamical equilibrium). The 
central BH was allowed to grow in mass by disrupting stars 
at R t .d. = (10~ 9 - 10~ 7 )^nb and fully accreting their mass. As 
Figure [6] clearly indicates, we can reproduce the evolution of 
such systems in a satisfactory manner using ME(SSY)**2. 
We have also checked that our results are insensitive to the 
particle number (as long as it is large enough) by repeating a 
few models with N p = 5 /V* instead of N p = N*. On the other 
hand, we have found the MC results to be more sensitive 
on the time step parameter than one might hope. For these 
sing le-mass models, fst = 0.005-0.01 gives the best results 
(see lFreitag & Ben3l2001l for an explanation of how the time 
steps are determined in ME(SS Y)**2; in rough terms fg t is a 
prescribed upper bound on St(R)/t r i x (R)). With larger values, 
the deflection angles in "super-encounters" become too large, 
leading to too little relaxation (and hence evolution) per unit 
of simulated physical time. 

The next step is to consider 2-component models in which 
a small fraction /heavy of the stars are significantly more 
massive than the rest with q = m ne avy /mught > 10. How- 
ever, there are no published results of this type using N- 
body simulations that can provide a well-controlled test case. 
For this reason we have undertaken our own A^-body sim- 
ulations using NBODY4, a code developed and made freely 
available by Sverre Aarseth 6 . Modifications were made to 
the code to include tidal disruptions and BH mergers. Over 
the years, Aarseth's Nbody family of codes have become 
central workhorses in a great numbe r of stellar dynamica l 
projects. They are described in detail in lAarsetlj ( fl999ll200l . 
NBODY4 can exploit a GRAPE board to accelerate t he com- 
putation by a very large factor jMakino et al.ll2003l) . which 
proved essential to obtain the results presented here. 

To represent stellar BHs in a population of age 5-10 Gyr 
q ~ 20-40 and /heavy — (1 - 3) x 10~ 3 would be adequate but 
such small fractions cannot be adopted usefully in A^-body 
simulations with A^ p < 1.3 x 10 5 , the highest number of par- 
ticles that can be used on the micro-GRAPE hardware at our 



6 http : / /www . ast . cam . ac . uk/~sverre /web /pages /nbody . htm 
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FIG. 7. — Comparison between our NBODY4 and ME(SSY)**2 simula- 
tions of a 2-component cluster model. The evolution of Lagrange radii for the 
indicated mass fractions of each component is plotted. The initial structure is 
an r\— model with r\ = 2 and fj, = 0. 1 . The population consists of light stars with 
flight = 1 Mg , shown with solid lines, and of 5 % (in number) of heavy ob- 
jects with mhcavy = IOMq, shown with dash-dotted lines. For these runs both 
types are treated as MS stars, subject to tidal disruption by the central BH. 
The physical scales are set by Rnb = 1 pc and N* = 64 000. The gray (orange) 
curves show the results of the W-body simulation, realized with N p = 64 000. 
The black curves are for the MC run which used Np = 640000. The N-body 
time unit is converted into FP unit assuming 7 C = 0.02. The Lagrange radii 
for the /V-body run are determined, at each snapshot, through a procedure 
of "orbital oversampling" in which the position of each particle on its orbit 
is sampled many times, with probability density dP/dR <x v~ (R), assuming 
a spherically symmetric potential centered on the IMBH. This way, one can 
follow a fractional mass as low as 0.001 which represents only 3.2 particles 
for the heavy stars. 



disposal. To have a reasonably large number of heavy par- 
ticles, we have chosen /heavy = 0.05 and q = 10 for a simu- 
lation with N p = 64000. The initial structure of the cluster 
is an 77-model with 77 = 2 and ^ = 0.1. For simplicity, we 
have assumed that all stars have a MS size and are tidally dis- 
rupted if they come within /? t .d. of the IMBH, itself treated as 
a massive particle (rather than an external potential). When a 
star is tidally disrupted its whole mass is given to the IMBH. 
The size is set to /?nb = lpc. MC models were run with 
N p = N« ("64k"), = 5N« ("320k") and N p = ION* ("640k") 
for higher resolution and to permit a better determination of 
the local density, particularly near the cluster center, as needed 
by the MC algorithm for robust results. Actually, the results 
turn out to depend very little on Ap. 7 . 

Figure Q offers a global view on how the spatial distri- 
butions of light and heavy particles evolve with time in the 
Af-body and MC simulations. For the A^-body simulation, 
the center of the system, from which distances are measured, 
was defined to be the (instantaneous) position of the IMBH. 
As the natural time scale is dynamical for the A^-body code 
(?nb) but relaxational for the MC algorithm (Tpp), one needs 

7 The N— body model was run at the Astronomisches Rechen-Institut in 
Heidelberg, on a PC equipped with a micro-GRAPE board. It required ap- 
proximately 2 weeks of computation. In contrast, 64k and 320k MC runs 
took about 0.5 and 4 hours on a 1.7 GHz laptop. 
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FIG. 8. — Comparison of the density profiles between our NBODY4 and 
ME(SSY)**2 simulations of a 2-component cluster model. Thick solid and 
dash-dotted lines show the mass density for light stars ("MS") and massive 
ones ("BH"), respectively; the thin lines are the total densities. The W-body 
and MC results are shown in gray and black, respectively. We also add 
straight lines representing power laws with 7 = 1.75 and 7 = 1.5. 



to specify j c to compare the results in the same time units. We 
find the best agreement with 7 C = 0.01-0.02, as was the case 
for the (DMBH-less multi-mass systems simulated bv Freitag 
et al. ( 2006). For the light particles, the concordance between 
the methods is excellent. The heavy particles, on the other 
hand, show some discrepancy. The MC code produces mass 
segregation at a rate almost equal to that seen in the N-body 
runs. The heavy objects appear to concentrate slightly more 
at the center before the whole cluster starts expanding slowly. 

The nature of the difference between the results from the 
two codes is seen more clearly in Figure|8]where a snapshot of 
the central density profiles at nearly the same time is shown. 
The MC run shows a Bahcall-Wolf cusp of BHs that extends 
all the way down to the resolution limit. In contrast, the 
A'-body profile appears to flatten slightly inside R w 0.01 7?nb- 
Given that the region with this flattened profile involves only 
5-10 BH particles at a time in the A'-body simulation, this 
mismatch could be deemed of little significance, if it were 
not consistently present in most snapshots. We have redone 
the MC simulations with or without large-angle scatterings or 
tidal disruptions of the MS stars and found that the results 
are not altered: in all cases, the BHs develop a slightly more 
pronounced innermost density peak than in the A'-body run. 
The fact that in the MC simulation the central BH is assumed 
to be fixed in position may be the cause of the difference; 
this is supported by the amplitude of the IMBH wandering in 
the A^-body run, ~ 0.01 7? NB (comparable to the spatial extent 
of the flattened profile). If this is the case, the effect should 
be less important in galactic nuclei, as far as the distribution 
of stars around the MBH is concerned because the wander- 
ing -essentially the manifestation of energy equipartition- de- 
creases with decreasing mass ratio m*/M. l lLaun & Merrittl 
120041 and references therein). In our A'-body simulation, this 
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FIG. 9. — Evolution of Lagrange radii for a "standard" MW nucleus model 
GN2 5. We plot the evolution of the radii of spheres that enclose the indicated 
fractions of the mass of various stellar species. Solid lines are for MS stars, 
short-dashed lines for white dwarfs, long-dashed lines for neutron stars and 
dash-dotted lines for stellar BHs. 



ratio is is ~ 10~ 4 and ~ 10~ 3 for light and heavy stars respec- 
tively. In a galactic nucleus with a 10 6 M Q MBH, the ratio is 
10~ 5 at most. 

Last we examine the density profiles shown in Figure [8] 
Specifically those obtained with the MC code (less affected 
by small-number effects) clearly indicate that the light ob- 
jects follow a profile compatible with 7=1.5 only at distances 
smaller than R\.s w 0.01 /?nb, whereas the radius encompass- 
ing a mass of stars equivalent to M. (an approximation to /?i n fl) 
is of order 0.2 -0.3 R^b- Only within R\.sis the number space 
density of stars dominated by BHs. Between R U 5 and R^n is 
a transition region in which 7 < 1 .5 for the light objects even 
though 7 ~ L75 for the heavy ones . Similar findings were 
obtained bv lBaumgardt et al.l(l2004bl) . Although these results 
do not invalidate the prediction from the Fokker-Planck treat- 
ment that light objects should form a cusp wi th 7 = 1.5 close 
to th e central (DMBH jBahcall & Wolflll977l Gnedin & Pri- 
mack T27)04ir ih"e'v indicate that, unless the fractional number 
of massive objects is unrealis ticallv high (as is the case in the 
test-computation presented by M urphy et al J 1 99 ll in their Fig- 
ure 9), this regime may only be attained in a very small cen- 
tral volume and therefore will be of relatively little relevance 
to real systems. 

4.2. Realistic Models 
4.2.1. Sgr A* -type models 

Next we consider models specifically intended to represent 
galactic nuclei. The parameters describing the initial condi- 
tions for these simulations are listed in Table|2] 

Let us first consider in some detail the evolution of our 
"standard model", run GN2 5 with 77 = 1 .5 and n = 0.05. These 
parameters are adapted to fit the observed enclosed mass 
profil e of the MW center (ISchodel et al.N2003HGhez et all 
120051) . The physics implemented in this model is limited to 2- 
body relaxation, tidal disruptions of stars by the central MBH 



(which accretes 50 % of the stellar mass) and direct plunges 
through the horizon. We use stellar population F. This is one 
of the highest resolution models with iV p = 8x 10 6 and each 
particle representing 26.65 stars (note that the MC code does 
not require a particle to stand for an integer number of stars). 

The overall evolution of the nucleus structure is depicted 
in thre e different (but essentially equivalent) ways in Figures 
l9l 1101 and 1111 In Figure |9] we present a general overview by 
showing how the Lagrange radii of the various stellar types 
evolve with time. The development of mass segregation is 
clearly apparent. Qualitatively, the region of influence of the 
MBH corresponds to the extent of the MS Lagrange radius 
for a fractional mass equal to the value of /1 = M./M c i, i.e. 
0.05. Deep in this region, the evolution is approximately ho- 
mologous. The stellar BHs concentrate in the center over a 
timescale ~ (1-2) X 10" 3 T FP w 4 - 8 Gyr. At the same time, 
the other stars slightly expand out of the center but the total 
density profile stays nearly constant. 

During this first phase, the BHs come to dominate the cen- 
tral mass density by forming a cusp around the MBH. This 
can be seen in Figure[H)] We note that, at late times, the cusp 
exponent becomes compatible with 7 = 1.75, but the lighter 
objects form a profile^itli^^l^. flatter than the Bahcall 
& Wolf i 19771) exponent. However, it must be stressed that, 
for this model, the stellar BHs never contribute more than 
~ 15 % of the number density in any region. Therefore they 
do not become a strictly dominant species, in the sense that 
they still experience most of their interactions with lighter 
objects. This is different from the situation studied bv Bah- 
call & Wolf (119771) . who only considered larger values for the 
number fraction /heavy of massive objects (their smaller value 
being /heavy =1/16 while we have /heavy = /bh — 0.002) and 
smaller mass ratios q = mheavy/'wiight (they have q < 10 while 
we have q w 20-30). Because our particle number is not 
large enough to treat the system on a star-by-star basis, it is 
still possible that, in a real MW-like nucleus, there would be a 
region very close to the central MBH in which the stellar BHs 
are numerically dominant and a clean IB ahcall & Wolfl(ll977l) 
cusp could form. Our results strongly suggest that the radius 
of this region is at least 100 times smaller than R[ n f\. 

All other stellar species react to the segregation of the stel- 
lar BHs by expanding away from the center. This evolution 
is very similar for all objects of mass significantly lower than 
that of the BHs, with the NSs showing slightly less expansion 
than the MS and WD stars. However, to the resolution limit 
of our simulations, the density profiles show no conspicuous 
central d epletion, such as a flatten ing or even a dip (as sug- 
gested bv lChaname & Gouldl2002l for pulsars around Sgr A*). 
Such a density decrease is apparent only in comparison with 
the initial conditions. It is very unlikely that this density de- 
crease can be revealed by observations in the Galactic center 
as a tell-tale indication of the presence of a cusp of stellar 
BHs. Also, MS stars of different mass es react essentially the 
same way, as can be seen in Figure^] and end up having the 
same density profiles. 

The fact that the stellar BH population is the main driv- 
ing cause for the evolution of the central parts of our nucleus 
models become s cle ar by running a simulation without any 
BHs (see Figure \12\ . The most obvious difference is that the 
overall evolution, now driven by the mass segregation of NSs, 
is of order 5-10 times slower, reflecting a correspondingly 
longer dynamical-friction time scale. The NSs are fully seg- 
regated only after of order 30^40 Gyr. Consequently, even a 
clear-cut observation that old visible stars form a 7 < 1 .5 at 
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Radius [pc] Radius [pc] Radius [pc] 

FIG. 1 1. — Evolution of the profiles of enclosed mass for GN2 5. The solid lines are the results of the MC simulation. For reference, he dashed lines show 
rj = 1.5 profiles adjusted on the total mass and half-mass radius of each component. The top thin line is the total mass, including the central MBH; it is compared 
to the observational constrains for the MW center (see Figurel2l. 



the center of the MW could not be interpreted as (indirect) 
evidence for the existence of a population of invisible BHs 
following a steeper profile: if BHs are not present, the system 
evolves too slowly to reach a relaxed state over 5-10 Gyr and 
the observed distribution may still reflect some "initial condi- 
tions" impose, for example, by a merger with another nucleus 
or by a large starburst due to massive gas inflow. 

We note that the choice of rj = 1.5 as initial condition 
is rather arbitrary. It is mostly motivated by the observa- 
tional constrains on the present-day stellar distribution around 
Sgr A* . We have considered models with rj in the range 1 .2 
(7 = 1.8) to 2.25 (7 = 0.75) to assess the importance of the 
initial density profile on the late-time structure and evolution 
of our models. In Figure^I we compare the evolution of two 
models that share the same physics and most initial condi- 
tions, including the total mass, the mass of the MBH, the stel- 
lar population and (approximately) the enclosed stellar mass 
within 1 pc, but different central profiles, namely rj = 2.25, 
corresponding to a shallow cusp, and our usual 77 = 1.5. The 



77 = 2.25 model shows more evolution in the first 6-8 Gyr as it 
"catches up" with the 77 = 1 .5 case. After t ~ 8 Gyr, however, 
both nuclei have similar structures. In both cases, the BHs 
dominate the mass density inside R ~ 0.3pc (where their den- 
sity is ~ 2 x 1O 5 M pc" 3 ) at t ~ 10 Gyr. At that time, the BHs 
and MS stars form cusps with 7 = 1 .7 - 1 .8 and 7 = 1 .3 - 1 .4, 
respectively (for R < 0.15pc) for both simulations. In other 
terms, in the region of influence of the MBH, a period of 
time of order f r i x (^infl) (which translates into 7- 10 Gyr for 
our MW-like models) seems enough to erase the details of 
the "initial conditions". The initial conditions of model GN2 5 
were chosen to be compatible with the overall mass distribu- 
tion in the Sgr A* cluster, as constrained by observations (see 
Figure^. In Figure^Jwe see that, despite mass segregation 
and the slight expansion of the lighter stars, the enclosed mass 
profile is still an acceptable fit to the Sgr A* data after 10 Gyr 
of evolution. This is primarily because the evolution amounts 
mostly to redistributing the various stellar types while keep- 
ing the total density nearly constant. It is evident that the ob- 
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FIG. 12. — Evolution of Lagrange radii for a nucleus model without stellar 
BHs (GN4 9). 
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FIG. 13. — Differential mass segregation amongst MS stars. For model 
GN4 4, in which the IMF extends down to 0.01 Mq (instead of being trun- 
cated at 0.08 Mq), we plot the evolution of Lagrange radii for MS stars in 
three different mass bins. The lightest objects expand only slightly more than 
the most massive ones. For each bin, we indicate the average mass (m*) and 
mass fraction f ra . 



servations of the current mass profile do not provide a strong 
constraint on initial nucleus properties, as long as they match 
the stellar mass enclosed within R ~ 2 - 3 pc. For the chosen 
initial conditions, the overall expansion of the cluster occurs 
on a time scale longer than the Hubble time but, as we will 
see, smaller nuclei expand significantly over a few Gyrs, ow- 
ing to their sho rter relaxation time (see Figure |2"U1) . 

In Figure fT31 we present the number of stars of various types 
within distances of 0.01, 0.1 and 1 pc of the MBH, as a func- 
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FIG. 14. — Comparison of the evolution of two nucleus models with similar 
initial conditions but different values of rj: r\ — 2.25 (GN20; top panel) and 
7] = 1.5 (GN17; bottom panel). Note that the conversion factors between 
/V-body and FP units (for length and time, respectively) and physical ones 
are different in each case. 



tion of time. It is again evident that it takes 3 -5 Gyr for the 
stellar BHs to concentrate in the inner pc. For a variety of 
rj values and stellar populations, we find that between 20 000 
and 30 000 of them populate this region after 5 Gyr. Without 
mass segregation their number would be of order 4-5 times 
lo wer. These numbers bracket the est imate of 25 000 obtained 
bv lMiralda-Escude & GouldfeOOCj). Sim i larly for a stellar 
population similar to our case S. iMorrisI ( fl99l found that 
some 3.6 x 10 4 BHs would dominate the stellar mass density 
in the inner 0.8 pc (see line 4 of his Table 1). This agree- 
ment could be taken as proof that the dynamical friction for- 
mula, used by IMorrisI ( 119931) and iMiralda-Escude & Gould 
J2000I) . captures the process of mass segregation quite accu- 
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FIG. 15. — Number of stars of different types within 0.1, 0.3 and lpc 
from the center for model GN2 5. Solid lines indicate MS stars, short dashes 
WDs, long dashes NSs and dash-dotted lines stellar BHs. The dotted lines are 
the result of the application of the analytical formula for dynamical friction 
(equationsl8landl9l for the BHs, assuming a static background defined by the 
jj-model and average stellar mass of the initial conditions. BHs that would 
have reached R = by dynamical friction are considered accreted by the MBH 
and not counted. 



rately. However we think that this agreement is actually rather 
fortuitous. In Figure^Jwe plot the predictions of the dynam- 
ical friction formalism, assuming circular orbits and a static 
background corresponding to the initial stellar distribution. 
BHs that reach R = are assumed to merge with the MBH 
and are not counted. Applied to our initial nucleus model, 
this computation overestimates the speed and magnitude of 
mass segregation. In particular, it leads to too many BHs be- 
ing accreted by the MBH and, consequently they lead to a 
fast decline in the number of BHs populating in the central 
region after t ~ 2 Gyr. For instance, from this simple treat- 
ment, one would expect only of order 7000 of them to inhabit 
the inner pc at t = 10 Gyr. As expected, this formalism also 
fails to reproduce the structure of the central BH concentra- 
tion by allowing BHs to sink in all the way down to R ~ 
and not taking into account their mutual interactions. Clearly, 
once the BHs dominate the mass density in some region, they 
start exchanging energy with each other at an important rate, 
a process which cannot lead to an overall contraction. Finally, 
based on the simple dynamical friction argument, one would 
erroneously expect all stars significantly more massive than 
the average, including the NSs, to segregate to smaller radii; 
this is clearly not seen in the numerical simulations. We note 
that using the local, self-consistent velocity distribution for an 
77-model instead of relying on a Maxwellian approximation to 
compute the dynamical friction coefficient makes a negligible 
difference. 

So far we have focused on our standard Sgr A* model. Ex- 
cept for mass segregation, its initial conditions were set to 
reflect the state of the MW nucleus at the present epoch in 
the sense that the enclosed mass profile (interior to ~ 3 - 5 pc) 
matches the observational constrains and that the stellar pop- 
ulation has a uniform age of 10 Gyr. Further stellar evolu- 
tion was not included. Such a model, chosen for its sim- 
plicity, is obviously not very realistic, not even entirely self- 
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FIG. 16. — Evolution of a model with stellar evolution and growth of the 
central MBH from a seed (GN7 8). The top panel shows the total number 
of stellar remnants in the nucleus. The total number of stars is 2.13 X 10 s . 
In the bottom panel, we show the evolution of the Lagrange radii for the 
various stellar species. A Lagrange sphere is specified by a fraction of the 
instantaneous total mass of stars of the corresponding species. The MBH 
grows from a seed of ~ 1000M Q to 3.95 X 10 6 Mq at t = lOGyr. Most of 
this increase comes from the accretion of a fraction 0.0653 of the gas emitted 
by stellar evolution. 
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FIG. 17. — Number of stars of different types within 0.1, 0.3 and 1 pc from 
the center for model GN7 8. 
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consistent. In particular, during the ~ 10 Gyr over which 
we allow dynamical evolution to proceed, the evolution of 
stars with ZAMS mass above ~ 1 M Q should be accounted 
for in principle. Also, the MBH may have significantly grown 
during such a long period. By considering a very different 
model, GN7 8, tuned to yield a Sgr A* -like enclosed mass 
profile after 10 Gyr of evolution, we show that the conclu- 
sions about the mass-segregation (and rates of interactions 
between stellar objects and the MBH, see next subsection) 
are largely insensitive to our assumptions about the past his- 
tory of the nucleus, within the framework of the assump- 
tions common to both models (spherically symmetry, evolu- 
tion in isolation, etc.). This conclusion also applies to the 
other, less radical, variations of the Sgr A* model that we 
have considered but do not discuss in detail. To help iden- 
tify models that are applicable to the Sgr A* cluster, in Ta- 
ble[3] we indicate the enclosed stellar mass M enc i(R) for R = 1 
and 3pc after 5 and 10 Gyr of evolution. Observations in- 
dicate that, for Sgr A*, M end (lpc) ~ 0.5-1 x 10 6 M Q and 
M end (3pc) ~ 0.5 - 1 x 10 7 M Q (see Fig.EJl. We note, that, be- 
cause the density of 77-models decreases steeply for R> Rb 
and we cannot afford large values of Rb, lest the central res- 
olution become insufficient, it is difficult to put enough mass 
within 3 pc of the center. 

Model GN7 8 is started as a cluster with t] = 3, i.e., no ini- 
tial central density cusp, containing a "seed" BH at its center, 
M.(0) = 10~ 5 M c i(0) (because 7 = but the velocity dispersion 
is isotropic, the few particles initially in the influence region 
are not in exact dynamical equilibrium). All stars are on the 
ZAMS at t = 0; as the simulation proceeds, they are turned 
into remnants at the end of their MS lifetime, according to 
prescription F of Tabled Natal kicks are imparted to NSs and 
stellar BHs. We setiV* =2.13x 10 8 , M cl (0) = 1.24 x 10 8 M Q 
and make the ad hoc assumption that 6.53 % of the gas emit- 
ted by stellar evolution is instantaneously accreted my the 
MBH, in order to get, at t = 10 Gyr, M. ~ 3.5 x 10 6 M Q and 
M./M c i ~ 0.05, similar to the parameters of most other mod- 
els. As tidal disruptions and coalescence also contribute to 
the growth of the MBH, we obtain M. = 3.95 x 10 6 M© at 
t = 10 Gyr. Because the central parts of the cluster strongly 
contract in the initial phase (see below), we had to simu- 
late clusters with different initial sizes to find a value that 
yield a good fit to the observed enclosed mass profile, namely 
tf NB (0)=16.2pc. 

We show the evoluti on o f the structure of this model in 
Fig. I16land plot in Fig. ^]the number of stars in the vicin- 
ity of the MBH. Nearly 90 % of neutron stars receive natal 
kicks strong enough to escape from the nucleus. A strong 
and relatively fast contraction of the inner regions starts at 
t ~ 3 Myr, which goes on, although at a much reduced rate un- 
til t ~ 100 Myr. This reflects the adiabatic contraction of the 
stellar orbits, nearly unaffected by relaxation on such a short 
timescale, in response to the growth of the MBH as it accretes 
the gas shed by massive stars turning into BHs and NSs. At 
t ~ 10 Myr the MS stars have formed inside R = 0.1 pc a pro- 
file compatible with the cusp p oc R ~ 2 predicted by the th eory 
for an initial distribution with 7^^3jOjmdan^UiljH?95^Fre- 
itag & Benz 2QQ2bfK The later evolution of the nucleus is 
again dominated by relaxation. The system of BHs reaches 
its highest concentration after some 2 Gyr. After that time the 
structure and evolution are essentially the same as that of the 
standard model. 

Our assumption about the fraction of the mass lost in stel- 
lar winds accreted by the MBH is ad-hoc. At early times 



it leads to highly super-Eddington growth (see bottom panel 
of Fig. II (St . It would be more physical to assume that the 
gas accumulates at the center until the Eddington-fed MBH 
can accommodate it but this would only introduce negligible 
changes in the results as long as this centra l gas reservoir is 
seen a s a point mass by the stellar system iFreitag & Benzl 
I2002bl) . In any case, the fate o f the interstellar gas i n a galac- 
tic nucleus is a complex i ssue (iDavid et al!l987alb[ Coker & 
Meria ll997lll999tlWilliams et al.lll999tlCuadra et afEool . 
well beyond the scope of this study. Because the MBH ac- 
quires the bulk of its mass on a timescale much shorter than 
the relaxation time but significantly longer than the orbital 
time of the stars affected by its growth, the results of our 
model apply to any situation of MBH growth respecting this 
hierarchy of tim escales, such as gas i nfall triggered by a galac- 
tic merger (e.g.. lBarnes & Hernauistll99 11119961) . 

To conclude the presentation of Sgr A* -type models, we 
take a look at the results for the rates of disruptive events. 
In Fig. we plot the accretion rates onto the MBH, for 
our standard, high-resolution model (GN2 5), for two lower- 
resolution models that include stellar collisions (GN4 5 and 
GN4 6) and for our "alternative" Sgr A* run with stellar evo- 
lution and progressive MBH growth (GN7 8). We plot the con- 
tributions of tidal disruptions (half of the mass of the star is 
accreted), coalescences, collisions (100 % of the gas liberated 
is accreted) and stellar evolution (for GN7 8). The mass lost 
in colli sions between MS stars is determined fro m our SPH 
results (iFreitag & Bendl2005tlFreitag et alJl2006t see §IXTV 
As for collisions between a MS star and a compact remnant, 
we considered two extreme assumptions: either we neglected 
them altogether (but counted their number), as in GN4 5, or the 
MS star was considered to be entirely destroyed in the process 
and the CR was left unaffected, as in GN4 6 (in another run, 
we assumed half of the mass of the MS star was accreted onto 
the CR). 

In Fig. [HJ] we show the number rate of tidal disruptions, 
coalescences and collisions for the same simulations. Irre- 
spective of the details of the models, we find that around 
t = lOGyr, the tidal disruption rate is \dN/dt\^ ~ 3-4 x 
10~ 5 yr _1 , in good agree ment with previous estimates for nu- 
clei of similar structure (j Rees 19881 iMagorrian & Tremainel 
11999b ISver & Ulmerll 19991) . Coalescences are less frequent 
by some 20-30 % but dominate the mass accretion rate owing 
to the important contribution of stellar BHs. In the models 
without stellar evolution, mass segregation is responsible for 
the significant increase in the coalescence rate taking place 
between t ~ 2 Gyr and t ~ 6 Gyr. The decline in the rates at 
later times is the consequence of the overall expansion of the 
nucleus. The contribution of collisions to the growth of the 
MBH never exceeds 10" 6 M Q yr _1 . At t = lOGyr, it is around 
5 x 10~ 7 Moyr _1 if CR-MS collisions are neglected and some 
4 times higher if these events are disruptive. As we will see in 
§ 14.2.31 collisions have also only little influence on the struc- 
ture of the galactic nucleus, as far as it can be resolved by our 
simulations. 

4.2.2. Models for Nuclei of Different Masses 

We now explore how nuclei less or more massive than our 
standard Sgr A* case evolve. We recall that, for a given rj 
value, we restrict ourselves to a one-parameter family of mod- 
els by keeping \i = M m jM c \ fixed and imposing R NB oc y/M„ 
a scaling is inspired by the M-a relation (see § I3.2> . Hence, 
the model is specified by t] and M.(0). We have considered 
M.(0) values ranging from 10 4 to 10 7 M Q and two values of 
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FIG. 18. — Rate of accretion onto the MBH from various channels. The top 
panel is for our standard model GN25. The middle panel is for simulations 
that include stellar collisions, GN4 5 and GN4 6. The lower (thick) curve for 
the contribution of collisions corresponds to run GN4 5 in which the collisions 
between MS stars and compact objects were neglected; the upper (thin) curve 
is for run GN 4 6 during which we assumed that such collisions always lead to 
complete destruction of the MS star. 100 % of the gas emitted in collisions is 
accreted by the MBH. The lower panel is for simulation GN7 8 which stalled 
with a seed BH (M. ~ 1000M Q ) and stars on the ZAMS. Simple stellar 
evolution is included with a fixed fraction of the gas emitted when a MS star 
turns into a remnant being accreted by the MBH (see text). For this run, the 
Eddington accretion rate (with 10 % conversion factor) is also plotted but the 
growth of the MBH was not limited by it. 
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FIG. 19. — Event rates for the same simulations as in Fig.EH 



rj: 1.5 and 2.0. 

We show the evolution of the model with M.(0) = 10 5 M Q 
and rj = 1.5 (model GN55) in Fig. [20] The most obvious 
feature is a faster evolution, when measured in years, than 
the Sgr A* nucleus, which simply reflects a shorter relax- 
ation time. This results in significant expansion of the clus- 
ter central parts over a Hubble time. For instance, the half- 
mass radius which showed hardly any change in the Sgr A* 
case, expands by a factor of about 2. In the models with 
M.(0) = 10 4 M Q , the whole cluster is expanding at t = 10 Gyr, 
with the half-mass radius of the rj = 1 .5 model increasing form 
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FIG. 20. — Evolution of Lagrange radii for a small galactic nucleus. The 
initial conditions for this model (GN55) are r\ = 1.5, [i = 0.05, a stellar popu- 
lation of type F, M, = 10 5 M Q andft NB =4.73pc. 



1.3 to 15.3 pc. The expansion proceeds like R(t) oc t" with 
j3 ~ 2/3 at large radii, as predicted for the s elf-similar ex- 
pansion powered by a central energy source llHenonl 1 1965t 
LShapiroll977tlMcMillan et al!l98lHGoodmanll984 Amaro- 
Seoane et al. 120041) . The central parts appear to expand slower 
with (3 ~ 1/2, a relation not yet explicitly reported in the 
literature, to the best of our knowledge, but consistent with 
the results of recent single-mass simu lations with the gas- 
dynam ical model of cluster dynamics jAmaro-Seoane et al.l 
12004 and NBODY4 (iBaumgardt et al.ll2004al) . Our exper- 
imentations with the gas model indicate that this is a loss- 
cone effect and that the transition from (3 ~ 1 /2 to 2/3 occurs 
around the critical radius. When tidal disruptions through 
loss-cone diffusion are prevented (i.e., stars are destroyed 
only when their energy becomes smaller than -GM,/R t ±), 
the whole cluster expands like R(t) oc f 2 / 3 . 

A straightforward consequence of this strong expansion of 
small nuclei is that one would have to start with initial condi- 

1 12 

tions much more compact to recover our assumed /?nb oc MJ 
relation between different nuclei in the present-day universe. 
However, this relation results from a naive application of the 
M-a relation. Observationally, a is a luminosity- weighted 
value integrated over the light-of-sight and averaged over an 
aperture covering a region much larger than t he one dynami- 
cally influenced by the MBH or by relaxation (iTremaine et al.l 
120021) . In fact, in no other case than the Milky Way is the 
region affected by r elaxation actually resolved by observa- 
tions (Merritt 2005). Furthermore, the M-a relation for 
M, < 10 7 M. is poorly constrained anyway. In future stud- 
ies of the evolution of small galactic nuclei, a larger variety of 
models should considered by allowing initial stellar clusters 
more and less dense than in our single-parameter families. 

In Fig.|^we make a direct comparison between the mod- 
els of different masses. We plot the evolution of the 0.3 and 
10% Lagrange radii for MS stars and stellars BHs and ac- 
cretion rates onto the central MBH (through tidal disruptions 



and coalescences). As we have seen, the distributions of WDs 
and NSs evolve similarly as that of the MS stars; they are only 
slightly more concentrated towards the center. For these plots, 
we have used "natural" units to stress the similarities between 
the various models. Time is expressed in Tpp, radius in 7?nb 
and mass in M c \(0). From the Lagrange radii evolution, one 
sees that models of different masses show a very similar struc- 
ture at the same value of f/7pp, which reflects the fact that 
evolution is driven by relaxation. Significant differences are 
only visible at small radii. They are the consequences of the 
"central boundary conditions" imposed by tidal disruptions 
and coalescences. Unlike relaxation, these processes intro- 
duce physical length scales in the system: R t ± and R$. The 
structure can only be independent of the size and mass of the 
model at distances larger than the corresponding critical radii. 

For the rj = 1 .5 series, the 0.3 % radius of the BHs contracts 
slightly faster at early times for more massive, larger nuclei 
with M, < 3.5 x 10 6 M Q . This seems to be the consequence 
of a bigger growth of the central MBH in the early evolution 
phase during which the stellar BHs segregate to the center 
(f < (1 - 3) x 10~ 3 7pp). In natural units, when the mass of the 
system is increased, the dynamical time at a given radius de- 
creases. For a fixed aperture of the loss cone, this would yield 
a higher accretion rate in the full-loss-cone regime (at large 
radii) and a larger critical radius while the empty loss-cone 
rate is unchanged. In fact, as we use a R oc M 1 / 2 scaling, the 
size of the loss cone, at a fixed R/Rnb value, varies approx- 
imately like #lc °c 

M. 1/6 for tidal disruptions (R ti . oc M 1 . 13 ) 

9 1 /2 

and like 0^ c oc MJ for coalescences (/? p i un ge oc M.). All this 
indicates that the coalescence rate should increase with M. in 
our families of models, as indeed is the case. The situation 
for tidal disruptions is more complicated, also because, as M, 
increases, a larger and larger fraction of MS stars are compact 
enough to withstand tidal forces down to the last stable orbit 
around the MBH. For instance, with M. = 1O 7 M , only MS 
stars more massive than w 0.6 M Q can be tidally disrupted on 
non-plunge orbits, i.e., have R t .d. > ^plunge- 

The "segregation phase" ends earlier and the concentration 
of stellar BHs is less pronounced in the massive nuclei. This is 
likely a result of the larger critical radius which yields an ap- 
proximately equal energy production rate (in "natural" units 
such as TV-body energy unit per Tpp) for a lower stellar den- 
sity in the inner regions. At the same time a larger critical 
radius explains the larger accretion rate, as stars are absorbed 
by the MBH when they are on wider orbits. To drift from 
large distances to these orbits the accreted stars have to dis- 
sipate less orbital energy and contribute less heating towards 
the stellar system. In Fig.|22j we verify that the energy pro- 
duction rate through tidal disruptions and coalescences with 
the MBH is indeed nearly the same for the different models 
with 77 = 2.0 d uring the expansion phase. As first realized by 
lHenoHfr975T) . during the gravothermal expansion of a cluster 
the conditions in the central regions where energy is produced 
are indirectly controlled by the large-scale structure. The lat- 
ter determines how much energy is transported outwards by 
2-body relaxation to drive the expansion; this "luminosity" 
has to be balanced by central energy production through dis- 
ruptions and coalescences, in way similar to hardening of bi- 
naries w hich powers the post-collapse expansion of globular 
clusters lShapiroll977tlInagaki & Lvnden-BelI198l: lHe gg iel 
U984tlHeggie & Hul2003l amongst others). 

4.2.3. Role of large-angle scatterings and collisions 
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FIG. 21. — Evolution of galactic nuclei of different sizes. The left panels are for models with r\ = 1.5, the right ones for 77 = 2.0. For all models, /1 = 
M.(0)/M c i(0) = 0.05. We consider MBH masses ranging from 10 4 and 10 7 Mq and scale the initial size of the cluster according to Eg. 1171 i.e., Rnb c< \/M m . 
In the top panels, we show the evolution of the Lagrange radii of fractional masses 0.1 and 0.003 for MS stars (solid lines) and stellar BHs (dot-dash). The 
triangles on the lower x-axis indicate 10 Gyr for the models in which this corresponds to less than 0. 1 in Fokker-Planck time units. In the bottom panels, we plot 
the accretion rate onto the MBH. Solid lines indicate the contribution of coalescences and the dot-dashed lines that of tidal disruptions (with 50 % of the mass of 
each disrupted star beeing accreted). 

The models withM. = 1O 4 M0 (3.5 X 1O 4 M0) is made up of only N p =N, = 6.1 X 10 5 (2.13 X 10 6 ) particles and the 0.003 Lagrange radius for BHs is determined 
with 3 (10) particles only, hence the large-amplitude oscillations. We used 4 X 10 6 (< N t ) particles for all other simulations plotted here. 



The discussion of direct collisions and large angle scatter- 
ings in galactic nuclei would warrant another, specific paper. 
Here we only consider the overall effects of these "strong 
encounters" between stars on the structure and evolution of 
galactic nuclei. 

In a few of our models, large-angle scatterings where in- 
troduced to test whether a significant number of stellar ob- 
jects are ejected from the cusp through this mechanism. In 
Fig. [23] we compare the number of ejections to that of coales- 
cences and tidal disruptions. For all stellar types, the number 
of ejections turn out to be smaller by a factor of a few. In par- 
ticular, stellar BHs are nearly 10 times more likely to me rge 
with the MBH than to be ejected. iLin & TremaineUl980T) ar- 
gued that it was some 30-70 % as likely for a cusp star to be 
ejected from the cluster as to be absorbed by the MBH. Our 
results do not confirm this estimate but the nuclei we consider 
are very different from the globular cluster model adopted by 



ILin & Tremaind in which the cusp is embedded in a large 
constant-density stellar core. These authors also argue that 
the probability of ejection from the core (as opposed to from 
the cluster) is 3-10 times that of absorption, so there is a 
possibility that large-angle scatterings would reduce the rate 
of coalescences and/or disruptions significantly even though 
they do not lead to numerous ejections. For our Sgr A* -type 
model, we find the number of BH-MBH coalescences to be 
reduced by some 35 -40 % by the effects of large-angle scat- 
terings. Other numbers of events are much less affected. We 
find no appreciable influence on the density profiles around 
the MBH at f w 10 Gyr. 

Collisions between stars are also found to have very little, 
if any, impact on our results concerning mass segregation and 
rates of coalescence and tidal disruptions. In most runs with 
collisions, we only considered collisions between MS stars. 
Encounters featuring one or two COs were counted but the 
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10" 3 0.01 

Time [FP units] 

FIG. 22. — Comparison between the rates of energy production due to tidal 
disruptions and coalescences with the MBH for nucleus mod els with r\ = 2.0. 
The data is for same models as in the right panels of Fig. 1211 The units for the 
energy production rate is the N-body energy unit per Fokker-Planck time, 
i.e., GM^j S^' B Tpp . The oscillations in the curves are essentially numerical 
noise. The different models show very similar energy production rates in the 
expansion phase. 



mass or trajectories of the stars were left unchanged. For a 
Sgr A* -type model, we find of order 6 x 10 4 MS-MS colli- 
sions in lOGyr (the number of collisions between particles in 
the simulation is smaller by a factor N^/Np = 53.3). The num- 
bers of MS-WD, MS-NS and MS-BH for the same period are 
about 3 x 10 5 , 2 x 10 5 and 2-4 x 10 4 respectively. Not sur- 
prisingly, collisions between compact stars are extremely rare 
and their number is therefore very uncertain given the reso- 
lution of our simulations. We find a number of WD-NS and 
WD-BH events of the order of 100-1000 (corresponding to 
only a few particle-particle encounters) and no collisions be- 
tween 2 compact stars. 

Despite the high velocities reached in the central regions 
(the median value of the relative velocity at "infinity" for 
colliding stars is about 500 kms -1 ), these encounters are not 
highly disruptive. Collisions very rarely result in mergers but 
the complete destruction of a MS star requires of order 20-30 
collisions if encounters with CRs are neglected; on average 
only about 4 % of the stellar mass is lost when two MS stars 
collide (see also lFreitag et alJl2004l) . We are not able to de- 
tect any significant effect of collisions on the central density 
profile of MS stars, down to a few 10~ 3 pc of the MBH, even 
when we assume that collisions with a CR result in the com- 
plete disruption of the MS star. This is somewhat surprising in 
view of estimates of the collision time such as done in Fig. [3] 
which suggest that inside ~ 0.01 pc of the MBH, a MS star 
should suffer from ~ 10 collisions in l OGyr. Strong central 
collisional depletion has been found bv lMurphv 
in their Fokker-Planck simulations but the nucleus models for 
which this effect was reported were much more collisional by 
construction than those used here. To determine whether col- 
lisions can have an observable effect on the stellar distribution 
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FIG. 23. — Cumulative number of events in simulation GN54. In the top 
panel, we show the numbers of tidal disruptions of MS stars and coalescences 
(i.e., plunges through the horizon of the MBH) for all stellar species. On the 
bottom panel, we plot the numbers of ejections from the nucleus, due to large- 
angle scatterings. Note that the number of ejections is always significantly 
smaller than that of coalescences with the MBH. 



at the Galactic center, simulations with a much higher resolu- 
tion are probably required. 

5. SUMMARY AND OUTLOOK 
5.1. Summary of simulations and results 

In this paper we report on our stellar dynamical simulations 
of multi-mass models of galactic nuclei. Our main goal was 
to investigate how stellar objects of different masses distribute 
themselves around a central massive black hole (MBH), in 
response to relaxation, a process known as mass segregation. 

This work is based on the use of a Monte-Carlo code which 
allows to follow the secular evolution of spherical stellar clus- 
ter in dynamical equilibrium over Gyrs. We performed about 
90 different simulations, to investigate the effects of vari- 
ous physical ingredients, assumptions about their treatment, 
of the initial nucleus structure and to perform some limited 
parameter-space exploration. For most models, 4 x 10 6 parti- 
cles where used, requiring a few days of computing time on 
a single-CPU PC. A few cases were computed with 10 6 or 
8 x 10 6 particles to establish that our results are not strongly 
affected by the limitations in numerical resolution. The latter 
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number is close to the maximum number of particles that can 
fit in the memory of standard 32-bit Linux PC. 

All runs included the effects of the gravity of a central 
MBH, of the self gravity of the stars, of 2-body relaxation, 
treated in the Chandrasekhar (diffusive) approximation, and 
of the tidal disruption of MS stars at the Roche limit around 
the MBH as well as direct coalescence with the MBH for 
stars too compact to be tidally disrupted. In most cases, stel- 
lar evolution was not included explicitly; instead the stellar 
population consists, from the beginning of the simulation of 
a mixture of MS and compact remnants corresponding to a 
single star formation episode that took place lOGyr ago. In 
a few models, explicit stellar evolution was included with all 
stars starting on the MS and turning into compact remnants 
(CRs) at the end of their MS life time. For simplicity, giants 
were not considered because, as far as mass segregation is 
concerned, only the mass of the star matters and the evolution 
of the stellar distribution, being a relaxational process requires 
timescales much longer than the duration of the giant phase. 
Stellar collisions and large-angle gravitational deflections (not 
accounted for in the diffusive treatment of relaxation) were 
considered in a small number of models. We made no attempt 
to determine whether a given star-MBH coalescence would 
occur as a gradual inspiral or a direct plunge through the hori- 
zon of the MBH. This question, of central importance for fu- 
ture low-frequency gravitational wave detectors such as LISA 
or BBO, has to be addressed in future work with more appro- 
priate numerical methods. 

All our runs are started as 77-models. They have a cen- 
tral power-law density cusp, p oc R v ~ 3 and steeper "cut-off 
at large radii, p oc R~ A . In most cases, we used parameters 
(mass of the MBH, stellar density around it, etc.) correspond- 
ing to the stellar cluster around Sgr A* at the center of our 
Galaxy. We did not try to reproduce the very peculiar spa- 
tial and age distribution of the bright IR stars observed within 
1 pc of Sgr A*. In this work, we adopt the position that these 
stars, useful as they are as probes of the gravitational poten- 
tial, are not representative of the overall stellar population at 
the Galactic center, assumed to be much older and therefore 
amenable to our treatment. The usefulness of this assumption 
is that it defines a well-posed problem which consists an inter- 
esting limiting case. Clearly other situations have to be con- 
sidered in future studies. For instance, an exciting scenario, 
in sharp opposition to our simplifying assumptions, is that the 
Sgr A* cluster and its central MBH have been formed progres- 
sively by the infall of rich stellar clusters, some of them con- 
taining intermediate-mass black holes (Hansen & Milosavl- 
ievic l2003UKim et al.ll2004 Idurkan & Rasioll2005l) . In this 
picture, which attempts to explain the existence of the very 
young, unrelaxed stellar populations and assumes the present 
epoch is not a special one, a stellar cluster should inspiral into 
the Galactic center everv few Mvr (but see Nayakshin & Sun- 
vaev T200^ forTrguments opposing this view and suggesting 
the young massive stars have formed in situ in an accretion 
disk). Such infalls would build up a mixed-age stellar popula- 
tion and reshuffle the orbits of stars already present in the nu- 
cleus quite significantly, and therefore yield a different mass- 
segregation structure. 

Based on our "standard" Sgr A* models and a somewhat 
naive application of the M—a relation, we have considered 
two families of galactic nucleus models of different masses, 
for M. in the range 10 4 - 1O 7 M . One family has 77 = 2.0, 
the other r) =1.5. The interval in M m was chosen mostly to 
cover the values that should yield gravitational wave signals 



in LISA band when a compact star inspiral into the MBH. 
We embarked in the present study as a first step towards more 
robust determinations of the rate and characteristic of such 
extreme-mass ratio inspirals (EMRIs). This range of models 
also covers systems that are both large enough (in terms of the 
number of stars) to be amenable to treatment with the Monte- 
Carlo method and small enough for relaxational effects to play 
a significant role over some lOGyr. 

To ensure that the Monte-Carlo code, based as it is on a 
number of simplifying assumptions, yield correct results, we 
performed a number of comparisons with simulations per- 
formed with highly accurate (but much more computationally 
demanding) direct-summation NBODY4 code. In particular, 
we performed a new NBODY4 simulation of a two-component 
model with a central massive object using 64000 particles. 
On the GRAPE hardware at disposal not more than ~ 10 4 
particles can be used; hence it is not yet possible to simulate a 
system with a realistic mass function using direct A^-body but 
this 2-component toy model demonstrated for the first time 
in a direct fashion that the MC code treats mass segregation 
around a MBH very satisfactorily. 

For the Sgr A* models, our main results are the following. 
In all cases, the stellar BHs, being the most massive objects 
(with a fixed mass of 10M Q or a range of masses, depending 
on the model), segregate to the central regions. This segrega- 
tion takes about 5 Gyr to complete. The nucleus then enters a 
second evolutionary phase which is characterized by the over- 
all expansion of the central regions, powered by the accretion 
of stellar mass (of very negative energy) onto the MBH. Al- 
though all species participate in the expansion, mass segre- 
gation continues in a relative fashion, as the system of BHs 
expands slower than the other components. The structure of 
the nucleus at distances from the center larger than a few pc 
is left unaffected by relaxation over a Hubble time. 

BHs dominate the mass density within ~ 0.2 pc of the MBH 
but we do not find them to be more numerous than MS stars in 
any region we can resolve (down to a few mpc, at t = 10 Gyr). 
Estimating the exponent for the density cusp the BHs form, 
p oc R~ 7 , is difficult because of numerical noise, but, in most 
cases, 7 is compatible with the Bahcall-Wolf value 7 = 1.75. 
In contrast, the less massive objects, such as MS stars, form 
a cusp with 7 generally in the range 1.3—1.4 which is sig- 
nificantly lower than the value of 1.5 nredicted bv Bahcall & 
Wolf (119771) .' This is also found in the 2-component A^-body 
simulation. After 5-10Gyr of evolution, we find of order 
2-3 x 10 3 ,6-8 x 10 3 and2x 10 4 stellar BHs within 0.1,0.3 
and 1 pc of the center, respectively. About 10 4 BHs coalesce 
with the MBH during a Hubble time. Using the formalism of 
the dynamical friction for objects on circular orbits in a fixed 
stellar background is an easy alternative for estimating the 
concentration of massive objects in the central regions. How- 
ever, for stellar BHs, although this approach offers a qualita- 
tively correct picture, it overpredicts the effectiveness of mass 
segregation. In the case of a model with 77 = 1 .5 (whose relax- 
ation time does not increase towards the center), this yields 
too large a number of BHs accreted by the MBH and too few 
being present within the inner 1 pc after some 10 Gyr. 

All types of objects lighter than the BHs, including the NSs, 
are pushed away from the central regions. Using the observed 
distribution of these object to infer the presence of segregated 
BHs does not seem to be possible though because, in the ab- 
sence of BHs, it would take the NSs more than 10 Gyr to form 
a Bahcall-Cusp of their own, even they would not receive any 
natal kick. 
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These results are not significantly affected by stellar col- 
lisions, large-angle scatterings, the initial 77 value. We also 
considered three different prescriptions for the masses (and 
types) of compact remnants and found no strong variations 
in the simulation outcomes. Most interestingly, an alternative 
model in which stellar evolution was included and the cen- 
tral MBH was grown from an IMBH seed (by accretion of 
an ad hoc fraction of the mass lost by stars when they turn 
into COs) yields basically the same structure of mass segre- 
gation (and same rates of coalescences end tidal disruptions) 
at t ~ lOGyr. These findings suggest that our main results are 
not very sensitive to the special "initial conditions" used, as 
long as they are fine-tuned to produce at t = 10 Gyr a given 
MBH mass and stellar mass within ~ 1 pc of the MBH. How- 
ever, it would be instructive to consider a larger variety of 
models in future work, including some with extended period 
of stellar formation. Our present assumption of a single burst 
of stellar formation maximizes the number fraction of stellar 
BHs and the time available for mass segregation. 

When large-angle scatterings (not accounted for in the stan- 
dard diffusive treatment of relaxation) are explicitly included 
(essentially as a special case of collisions), they are found to 
have only little impact on the rate of tidal disruptions or coa- 
lescences with the MBH. A stellar BH is about 10 times more 
likely to by swallowed by the MBH than to be ejected from 
the nucleus. In contrast with this, in th eir multi-mass TV-body 
simulations, Baumg ardt et alJ ll2004rJ) find that all stellar BHs 
except one are ejected from the cluster and ascribe this result 
to strong interactions with objects (generally another stellar 
BH) deeply bound to the IMBH. These interactions are likely 
to be "resonant", i.e., the three objects (including the IMBH) 
form a strongly interacting, chaotic configuration for many 
orbital time until one of the lighter objec ts is ejecte d (Baum- 
gardt, personal communication; see e.g.. lHutl|l993l) . In prin- 
ciple, this mechanism can be included into the MC code by 
extending the loss-cone treatment used for tidal disruptions 
and coalescences to interactions with the binary consisting of 
the MBH and the most bound stellar object and resorting to 
explicit integration of 3-body motion when a close interaction 
between the binary and a third object is deemed to occur. 

However this process is probably of little importance in 
galactic nuclei as a rough analysis suggests. Let's write E e j 
and S p ] to denote the cross sections for strong interaction 
with the innermost stellar object (followed by ejection from 
the nucleus) and for direct plunge through the horizon of the 
MBH, respectively. Then assuming t hat the scaling laws for 
three-body interactions established bv He ggie et al.ldl996l) ap- 
ply all the way to mass ratios as extreme as considered here, 
we estimate E e j/E p i ps (a/Rs)(m/M,) where a is the semi- 
major axis of the stellar object deeply bound to the MBH, 
Rs, is the Schwarzschild radius of the MBH and m is typi- 
cal mass of stellar objects. Now, for the interaction to be 
resonant, the inner binary has to be well separated from the 
other objects in the cusp, a < where ^1* is the radius con- 
taining (on average) one stellar object. Assuming a power- 
law density cusp n oc R' 1 inside the influence radius R m ^ 
of the MBH, one finds R u ps R^m/M.) 1 /^ . Therefore 
if the stellar velocity dispersion outside is a, one finds 
Sej/Sp! ps (c/a) 2 (m/M.) ( ' H ' )/(3_7) . For a a = 20kms _1 glob- 
ular cluster containing an 1O 3 M IMBH with 7 = 1.5, this 
ratio is of order 10 5 . But this is reduced to 10~ 2 for a galactic 
nucleus with M, = 1O 6 M and a = 200 km s" 1 . 
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FIG. 2 4. — C omparison between the distribution of transient X-ray sources 
found by Muno et al. 1 2005b) at the Galactic center and the results of one of 
our simulations (GN25, at t = 9.19Gyr). The observational data are repre- 
sented by diamonds connected by dotted lines. The smooth curves, one for 
each stellar species are the simulation data. Plotted is the number of sources 
whose sky position projects within a given distance of the center of the nu- 
cleus. This number is normalized to 1 at Z? n orm = 23 pc. A distance to the 
Galactic center of 8 kpc has been assumed. The 7 seven transients are more 
concentrated around Sgr A* than any stellar component in the simulation. 
See text for an assessment of the statistical significance of this result. 



5.2. Astrophysical applications, including future work 

Although we have not attempted a realistic modeling of the 
Galactic center, it is tempting to apply our results to one spe- 
cific observation of the Sgr A * region. U s ing the Chandra 
X-ray space-born observatory, iMuno et all (l2005bl) have de- 
tected 7 transient sources which appear to be much more con- 
centrated around Sgr A* than the overall stellar population. 
Here we examine whether this may be a direct consequence 
of mass segregation, if these sources are all stellar BHs ac- 
creting from a lower-mass companion. We make the strong 
assumption that these binaries are not formed or affected by 
interactions with other stars such as 3-body binary formation, 
partner exchange, ionization, etc. Instead we consider they 
just react to 2-body relaxation as point objects with a total 
mass app rox imated by the mass of the stellar BH. 

In Fig. [24] we perform a graphical comparison between the 
observed distribution of X-ray transients and the distribution 
of the various species, most importantly the BHs, in our high- 
resolution simulation GN2 5 at t = 9.29 Gyr. Clearly, the tran- 
sients are more centrally concentrated than the BHs in the 
simulation but given the small number of observed sources the 
plot itself is not sufficient to rule out our naive model for their 
distribution. If we pick up at random 7 sources with projected 
distance from the center smaller or equal to 23 pc following 
our "theoretical" BH profile, we find that their distribution is 
at least as concentrated as the observed one (in the sense of 
the Kolmogorov-Smirnov test) in some 15 % of the cases. It 
is therefore at this point not possible to exclude that the tran- 
sients owe their peaked profile purely to mass segregation but 
this seems somewhat unlikely. 

As pointed out by Mu no et all (l2005bl) . the rate of binary 
interactions should also increase steeply towards the center 
and this probably combines with mass segregation to produce 
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the observed distribution. In comparison with the siti ation man & Alexander 



in a globular cluster with a well-defined core velocit ' dis- 
tribution, the problem of binary dynamics in the vicin ty of 
a MBH is complicated by the fact that there is no clear-cut 
definition of the hard-soft transition. The Keplerian veloc 
ity dispersion increases virtually without bound when 01 e ap 
proaches the center. This may affect a binary on an or lit of 
relatively large semimajor axis a around the MBH becai se 2 
body relaxation will cause the orbit to reach down to a 



value 

-?)) 
! dis- 
sper- 
e., at 



R perl = (1 - e)a «C a over a timescale of order f r i x ln( 1 /( 
(e.g.. lFrank & Reeslll976l) . Therefore, a binary may b 
rupted even if it is hard relative to the local velocity d 
sion at the position where it spends most of its time, 
distances of order a from the MBH. 

The most extreme type of dynamical interaction a rJinary 
can experience is the tidal separation of its members if ts or 
bit brings it within ~ abinC^»/ m bin) 1//3 of the MBH. Hei 2 
is the semimajor axis of the binary itself and nibin its nass. 
This process is of great interest by itself both as a way t > ere 
ate "hypervelocit y stars" and to deposit a star on a ti ght orbit 
around t he MBH JffiIIslll98i lYu & Tremainel2003t Gilalan 
dris et al 

Given a model of the stellar distribution, one could est mate 
an average lifetime for a binary of known properties and hemi- 
major axis (assumed fixed), accounting for the low-/? p L ex- 
cursions caused by relaxation. However the complex question 
of binary dynamics in a galactic center would be better treated 
through self-consistent stellar dynamical simulations »f the 
sort presented here but including binary processes. Henon- 
style Monte-Carlo code are particularly well suited for fol- 
lowing the evolution of large systems with a significant] frac- 
tion of binaries whose interaction can be computed accurately 
by direct 3- and 4-body integration, as has already bei n re- 
alized in the context of voung and globular clusters (Giersz & 
Spurzem l2003HFregeau et al.l2005tlGurkan et al.l2005l) 

Concerning the prediction of EMRI rates and prop* rties, 
the determination of how 2-body relaxation shapes the i tellar 
distribution around the MBH is only a first -crucial- st sp. A 
robust estimate of the fraction of stars that eventually nspi- 
ral into LISA band, rather than plunge directly throuj h the 

e the 
;ntric 
2inis- 
than 



horizon while still on a wide orbit will probably requi 
development of a specific code. For stars on very eco 
orbits, one needs to follow the combined effects of GW 
sion and relaxation on a timescale significantly shortei 
allowed by the present ME(SSY)**2 code where the time 
steps are a function of the distance from the center and cannot 
depend explicitly on orbital parameters, lest conservation of 

energy become impossible. 

The work of iHopman & Alexanderl l|2p05) (itself in; pired 
bv lHils & Benderll 19951) indicates a promising avenuel The 
vast majority of EMRIs enter the GW regime when theii orbit 
is confined deep inside the region of influence of the r 1BH. 
Hence one could develop a code specialized in the d) nam- 
ics of stars on quasi-Keplerian orbits around an MBH. iHop- 



consistently using 
the Henon Monte- 



have followed the secular change of ec- 



centricity and semimajor axis of individual stars assuming a 
fixed given stellar Dackground to determine the diffusion co- 
efficients for 2-boi ly relaxation. A powerful development of 
their method woulp be to evolve the stellar distribution self- 
a treatment of relaxation borrowed from 
Carlo approach. One would use individ- 
ual time steps to b;tter follow the evolution of stars on high- 
eccentricity orbits mtil their fate (inspiral, plunge or, possibly 
ejection) makes no more doubt. A Monte-Carlo code could 
easily cope with thg 10 6 - 10 7 stars within the influence radius 

on a star-by -star b; sis. 

Recenl l\ iHopm in & Alexanderl fcOOfil) have considered, 
for the first time in the study of EMRIs, the role of "reso- 
nant relaxation", i e., of the random changes in eccentricity 
and orientation of Ihe orbital planes due to the non-vanishing 
but fluctuating ton [ue exerted on an orbit by the other orbits, 
each c onsidered as an elliptical mass wire C Rauch & Tremaina 
These authors find that resonant relaxation can in- 
crease the EMRI rkte by of order a factor 10, an exciting re- 
sult which is calling for confirmation by other computation 
techniques. Unforjunately, although strictly also a 2-body ef- 
fect, resonant relaxation is unlikely to be amenable to the type 
of local 2-body inl eractions at the core of the Henon Monte- 
Carlo method. 
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n 9 1 30 

U.Z1 JZ 


i 

4 


0.041 


1L,L 




GN 6 4 


2.0 


0.05 


p 


1 




4 


1 1 23 


TC C 




GN 6 5 


2.0 


0.05 


p 
r 


J.J 


9 1 39 
Z. 1 JZ 


4 


9 1 n 
zi .u 


TP P 






2.0 


0.05 


p 


i n n 

1U.U 


A flQfl 
o.uyu 


A 
4 


IS so 

JJ.JU 


TP P 




GN 6 7 


1.5 


0.0283 


c,, 
r 11 


J.J 


9 1 39 
Z. 1 JZ 


A 
4 


98 n 

ZO.U 


tp wfNir\ 


/SE - U.UZ3 


GN 6 8 


1.5 


0.0566 


r LI 


j. J 


1 .uoo 


1 
1 


1 7 f*A 
1 / .04 


TP QTh 






1.5 


0.0566 


Fu 




1 066 


] 


17 64 


TP ^FfNK^ 




GN 7 


1.5 


0.05 


p 


J.J 


9 1 39 
Z. 1 JZ 


A 
4 


98 n 

Zo.U 


TP P 


JU /C IVlo ITldSS dCLlcLcU OIlLO Drl 111 IVlo-Orl LOlllSlUllS 


GN 7 4 


1.5 


0.0566 


R^n 
DOU 


j. J 


1 .WOO 


1 
1 


1 7 f^A 
1 / .04 


TP ■sF 




GN7 5 


1.5 


0.0566 


R^n 
Doll 


J.J 


1 .uoo 


1 
1 


1 7 f^A 
1 / .04 


TP 'sFriSJK'^ 




GN7 6 


3.0 


io- 5 


Fu 


r^y 


2.132 


4 




1 C, on(JNrL) 


, 6.53 % mass from SE accreted by MBH 


GN7 7 


1.5 


0.0566 


BSu 


3.5 


1.066 


1 


10.0 


1 C, irl(iNlv) 




GN7 8 


3.0 


io- 5 


ru 


r** tj 


9 1 19 

Z. I jZ 


l 

4 


1 £ 9 

10. z 




, 6.53 % mass Irom SE accreted by MBH 


GN7 9 


1.5 


0.0566 


BSu 


3.5 


1.066 


1 


7.0 


TC, SE(NK) 




GN80 


1.5 


0.1 


BS 


3.5 


1.066 


1 


17.64 


TC 




GN81 


1.5 


0.045 


BSu 


3.5 


1.066 


1 


8.0 


TC, SE(NK) 




GN82 


1.5 


0.045 


BSu 


3.5 


1.066 


1 


5.0 


TC, SE(NK) 




GN83 


1.5 


0.05 


F 


0.035 


0.02133 


2.133 


2.8 


TC,C 




GN84 


2.0 


0.05 


F 


0.035 


0.02133 


2.133 


2.1 


TC,C 




GN85 


1.5 


0.05 


F 


0.01 


0.006094 


0.6094 


1.50 


TC,C 




GN8 6 


2.0 


0.05 


F 


0.01 


0.006094 


0.6094 


1.12 


TC,C 




GN87 


1.5 


0.045 


BSu 


3.5 


1.066 


4 


8.0 


TC, SE(NK) 




GN88 


2.0 


0.05 


F 


0.035 


0.02133 


2.133 


0.5 


TC,C 




GN90 


1.5 


0.1 


F* 


3.5 


1.066 


4 


15.0 


TC,C, SE 


Initial age of stellar pop. 5 Gyr 



NOTE. — A^*(0) is the initial number of stars (in 10 8 ). N p is the number of particles (in 10 6 , generally A^ p <C N*). ^nb is the N- body length scale. If not indicated otherwise, 
the time step parameter is f§ t = 0.04 and the Coulomb logarithm is In A = In 7 € jV* with 7 C = 0.01 . 50 % of the mass of tidally disrupted stars is accreted onto the MBH. 
(a) p. "pi ( i uc j a r ,, ; BS: from Belczynski, solar metallicity; BP: from Belczynski, metal-poor (see text and TableFTl. An "u" indicates an unevolved IMF. An "*" indicates a pecularity 
explained in the "Comments" column. 

C: collisions; LA: large-angle scatterings; SE: stellar evolution (NK: natal kicks); TC: tidal disruptions and coalescences with MBH 
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TABLE 3 

Stellar Mass Enclosed by R = 1 pc andR = 3 pc at 5 and 10 Gyr 



— 

Name 


^ 

Time 

(Gyr) 


^Wi(lpc) 
(1O 6 M ) 


A/f 11 „ \ 

™encl(3pc) 

(10 6 M o ) 


— ^ 

Time 

(Gyr) 


«encl(lpc) 
(10 6 M Q ) 


«encl(3pc) 

(10 6 M Q ) 


GNO 4 


4 97 


2 82 


13 00 


10 00 


2.04 


1 1 00 


GNO 6 


5 03 


71 


4.40 


10 02 


0.31 


3 58 


GNO 7 


5.16 


86 


4.61 


10 03 


0.54 


4.04 


GNO 8 


5 01 


83 


4.55 


9 98 


0.45 


3.85 


G N 1 


4.94 


0.83 


4.56 


10.04 


0.41 


3.76 


G N 1 1 


5 02 


76 


4 50 


10 02 


0.42 


3.86 


GN 1 2 


5 01 


97 


5 59 


9.98 


0.46 


4.59 


GN1 3 


5.04 


0.73 


5.25 


9.99 


0.34 


4.31 


GN 1 4 


5 00 


85 


8 39 


10 05 


015 


6.93 


GN 1 5 


6.43 


1 21 


5.24 


7 67 


1 18 


5.15 


GN 1 7 


S OQ 

J.U7 


U.70 


H.J l 


Q 19. 
y. 1 o 


n 77 


j.oy 


GN 1 8 


5 04 


1.25 


5.45 


10 03 


95 


4.61 


GN 1 9 


5.04 


1.12 


5 80 


10 04 


87 


4.91 


GN2 


4.95 


1.16 


6.28 


10.08 


0.92 


5.33 


GN2 1 


4 91 


1 04 


4 65 


10 01 


0.82 


4.06 


GN2 2 


4 88 


2 70 


1 1 on 


Q Q4 
y.y^ 


2 11 


7. JO 


GN2 5 


5.05 


97 


4.49 


10 04 


75 


3 82 


GN2 6 


4.90 


2.05 


8.14 


9.97 


1.64 


7.14 


GN2 9 


5 OS 


99 


4 51 


9 89 


0.80 


3.91 


GN3 


5.09 


0.93 


3.86 


10.14 


0.71 


3.26 


GN3 3 


4 93 


1 05 


4.67 


10 21 


0.83 


4.07 


GN3 4 


4 96 


95 


3 89 


9.96 


0.74 


3.32 


GN3 6 


4.85 


0.87 


4.30 


9.99 


0.67 


3.57 


GN4 4 


5 01 


98 


4.54 


10 03 


76 


3.89 


GN4 5 


5.29 


1.04 


4.62 


10.03 


0.82 


4.06 


GN4 6 


5.1 1 


1.05 


4.65 


9.69 


0.85 


4.11 


GN4 8 


4.94 


0.72 


2.98 


10.1 1 


0.54 


2.43 


GN4 9 


6.43 


1.19 


5 00 


7 67 


1.16 


4.89 


IjlM u 


5 04 


92 


4 14 


10 11 


71 


3 37 


GN5 3 


4 87 


1 04 


4.67 


10 12 


81 


3.97 


GN5 4 


5 07 


98 


4.54 


9 90 


76 


3 87 


GN5 5 


4 98 


065 


38 


10 01 


041 


25 


GN5 6 


5.00 


0.19 


1.07 


10.07 


0.14 


0.78 


GN5 7 


5 00 


44 


2.24 


10 04 


33 


1 78 


GN5 8 


5 00 


1 05 


4 58 


9 98 


83 


3 95 




J.Uj 


1 64 


7 30 


10 58 


1 56 


O.oo 


GN 6 


4 98 


1 06 


4.67 


9 82 


93 


4.18 


GN62 


5.01 


0.079 


0.46 


9.98 


0.048 


0.29 


GN63 


5.01 


0.25 


1.39 


9.95 


0.17 


1.01 


GN 6 4 


4 98 


54 


2.95 


10 04 


0.40 


2.30 


GN 6 5 


4 96 


1.16 


5 85 


10 06 


0.96 


5.02 


GN 6 6 


5 14 


1 30 


8 06 


9 91 


1.38 


7.88 


GN 6 7 


4 93 


0.64 


2 22 


10 01 


0.51 


1.90 


GN 6 8 


5 06 


74 


2 81 


9 87 


57 


2.32 


GN 6 9 


4 95 


72 


2 74 


10 02 


0.55 


2.27 


GN7 




1 06 


4 68 


Q 7 1 
y. 1 1 


U.OJ 


4 10 


GN7 4 


5.09 


0.80 


3.01 


9.90 


0.63 


2.54 


GN7 5 


5 02 


77 


2 94 


10 16 


59 


2.44 


GN7 6 


4 87 


53 


2 23 


10 22 


44 




GN7 7 


5 01 


1 20 


4 83 


10 09 


91 


3.96 


GN7 8 


4.90 


0.99 


5.01 


9.99 


0.79 


4.14 


GN7 9 


4.93 


1.62 


6.51 


10.04 


1.17 


5.20 


GN80 


4.81 


1.01 


4.01 


10.07 


0.79 


3.48 


GN81 


4.99 


1.25 


5.49 


10.00 


0.94 


4.41 


GN82 


5.01 


3.69 


17.00 


10.05 


2.68 


13 


GN84 


5.01 


0.025 


0.15 


10.02 


0.015 


0.087 


GN85 


4.99 


0.0055 


0.033 


10.02 


0.0031 


0.018 


GN8 6 


5.00 


0.0061 


0.036 


10.03 


0.0032 


0.019 


GN87 


5.01 


1.23 


5.42 


10.01 


0.90 


4.31 


GN88 


5.00 


0.0075 


0.042 


9.98 


0.0037 


0.022 


GN90 


5.00 


2.36 


9.19 


10.00 


1.86 


7.91 



